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ABSTRACT 

The purpose of this contract was to determine the effect of the synthesis 
approach (single or multirate) on the machine requirements for a digital 
control system for the space shuttle boost vehicle,, The study encompassed 
four major work areas: synthesis approach trades, machine requirements 

trades, design analysis requirements and multirate adaptive control techniques. 

The primary resuits of this study are two multirate autopilot designs for the 
low Q and maximum Q flight conditions that exhibits equal or better performance 
than the analog and single rate system designs. Also, a preferred technique 
for analyzing and synthesizing multirate digital control systems is specified. 
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1.0 INTRODUCTION 

Single rate digital autopilot designs have been plagued by the divergent require- 
ments of sampling at a high rate to prevent folding unwanted signals into the band- 
width of the control system ^ile also trying to keep the sampling rate at a mini- 
mum to avoid finite accuracy arithmetic problems. The solution to this problem is 
to go to a multirate design where the sampling rates for the various computations can 
be specified in a more consistent manner. For example ^ cascade mulrirate Filters can 
be defined with a fast update rate front end to circumvent the folding problem and 
to allow optimal performance of the high frequency filter requirements. The low 
fiw|uency back end filter can be defined to alleviate arithmetic problems while 
meeting the low frequency filter requirements. Similarly, multiloop multirate de- 
signs can be used to meet different sampling rate requirements in different loops. 

1.1 Study Approoch 

This study was performed to determine muitirate analysis techniques and to study the 
effect of the system synthesis approach (single or multirate) on the machine require- 
ments for a space shuttle type boost vehicle at low Q and maximum Q flight con- 
ditions. The applicability of adoptive control techniques to space shuttle was also 
. investigated . 

The study nativities can be grouped into the following major work areas: 

Synthesis Approach Trades 
Machine Requirement Trades 
Design Analysis Requirements 
Multirate Adaptive Control Techniques 

This report is divided into sections covrring these four tasks. A summary of the 
major study activities is given below. 

The Synthesis Approach Trades task is irrvolved with applying synchronous and multi- 
rate sampled data analysis and realization techniques to two fixed point in time 
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linearized models of the plant and control system furnished by NASA-JSC. First, 
the models were verified on the Boeing analysis program MDELTA and then the 
analog control system was digitized using Tustin's bilinear transformation with fire- 
quency prewarping. Next, o digital control system was designed starting with a 
z-^nsformed model of the plant. Lastly, a multirate control system was designed. 
Comparisons were made between the gain and phase margin characteristics for the 
different designs. Hybrid simulations of the systems were implemented and the 
analysis results were verified using these simulations. 

The objective of the Machine Requirements Trades task was to determine the effect 
of the synthesis approach on the machine requirements. The primary machine re- 
quirements and the effect of sampling rate on these requirements are reviewed. The 
machine requirements for a synchronous implementation versus a multirate implemen- 
tation are discussed. 

The Design Anoiysis task involved studying techniques for analyzing multirate sampled 
data systems. The primary existing techniques are reviewed and the state variable 
approach is discussed in detail. This latter approach was used to perform the analy- 
sis work in Task I. 

A literah re search and study of various parameter adaptive control schemes were 
conducted. Extensions to these techniques to moke them applicable to a shuttle 
type application are outlined. The pseudo adaptive autopilot used on the SRAM 
missile is discussed. Two possible digital implementations for these adaptive control 
systems are discussed. 

1 ,2 Study Results 

The contract findings are summarized in the following paragraphs: 

Synthesis Approach Trodes - Single rate and multirate yaw plcme autopilots were 
designed for the two fixed point in time plant models. For the low Q system the 
single rate and multirats designs gave near identical performance when judged in 
terms of frequency response characteristics although the imiltirate design exhibited 
slightly better phase stabil izotion of the first bending mode. From a frequency 
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mponse standpoint. Hie max Q muiHrats design gave much better gain stabilization 
of Hie body bending modes than the single rate system. The hybrid simulation re- 
sults verified Hie analysis predictions In all cases considered. A multirate design 
9 methodol'igy was developed and is outlined below: 

1* Continuous root locus 

2. Z-plane root locus and frequency response at sampling rates 
of interest 

3. Determination of high update rate and low update rate filter 
requirements 

4. Multlrote frequency response and root locus 

5. Fine tuning multirate design based upon (4). 

Machine Requirements Trades - The effect of the single rote and multirate designs 
resulting from Task 1 on ‘ie machine requirements was negligible from a computer 
wordiength standpoint. This is primarily due to the fact that the allowable range 
in variation of the sampling period was restricted to a positive interger times T = .04 
sec. which Is Hie basic frame time for the Space Shuttle G&C computer. It was 
determined from both folding «id stability standpoints that the lowest acceptable 
sampling rate had a period of T = .08 seconds. The frame time of .04 sec. was 
slow enough to begin with to preclude any finite accuracy arithmetic problems for 
the filters considered for this system. The multirate design did free up computer 
time over the T = .04 single rate designs os not all the computations needed to be 
performed every cycle. Also, the max Q multirate design ./js a simpler design to 
implement because it was possible to replace a first order filter by a sample and 
hold element. The hybrid con^uter wordiength studies verified *’he aforementioned 
conclusions os they showed that all filters could be implemented on a 16 bit computer. 

Design Analysis Requirements - The two primory techniques that were studied for 
an-Jyzing multirate systems were Hie switch decomposition method of Krone [11] 
and the state transition approach of Kalman and Bertram [3]. The switch decompo- 
sition technique b shown to present severe computer implementation problems in the 
area of separating the branch varidsles from the loop equations. The state variable 
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approach is domonstrated to be straightforward to apply to a wide class of sampled 
data systems and to be easily implemented on a digital computer. The state transition 
method Is therefore recommended for analyzing multirate systems and a preferred 
method of defining the state transition matrices is specified. 

Multirate Adaptive Control Techniques - T'le parameter adaptive control algorithm of 
Lion [22] is shown to offer the best chance for doing online adqstive control based 
upon Identifying the vehicle characteristics. Pseudo-adaptive control techniques as 
used on the Boeing SRAM are also outlined and took as if they would offer lr:s 
technical risk than a full adaptive system. 
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2.0 . NTHESIS APPROACH TRADES 

The purpose of this section is to outline a methodology for designing multirate 
compensation for sample data control systems. To ?eet this goal, the section is 
broken into seven subsections. Subsection 2J presents two analog plant models 
which ore used as examples in developing the design methodology. Subsection 2.2 
addresses how different sampling rates c.. i affect the control system performance. In 
Subsection 2.3 synchronous compensation is formed by applying Tustin's method with 
frequency scale pre-warping to the analog compensation in the two plant models; 
Subsection 2.4 presents a single rate digital design for one of the plant models where 
the analog compensation has been removed and replaced with a sampler and zero-order 
hold. Subsection 2.5 presents a multirate design for both plant models. In Subsection 
2.6 the results of a simulation which verify the multirate designs are presenied. Sub- 
section 2,7 summarizes Section 2.0 with a general outline on the methodology tc use 
In designing multirate systems. 


I I 


2.1 Continuous Plont Gain ond Phase Margins 

The block diagrams of the two plant models used In developing the digital filters are 
shown in Figures 2.1.1 and 2.1.2. Both systems model the yaw dynamic of the space 
shuttle. The first system is for the lift off or low Q flight conditions. The second 
system is for the max Q flight conditions. In each ccse the design criterion is the 
some. The digital compensation is to replace the analog compensation while meeting 
the following constraints: 

1) The rigid mode gain and phase margins are to duplicate those of the 
continuous system. 

2) The first body bending mode is to be phase stabilized. 

3) The higher body bending modes are to be gain stabilized. 

The Nichol's chart for each of these systems where the loop has been ope r<ed at K is 
shown in Figures 2.1.3 and 2.1.4. 
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Figure 2.1.3: Continuous System Nichols Chart for Lift Off Flight Condition 
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2.2 Scwnpling Rate Requiremenh 

The choice of the sampling frequency is of prime importance when convertirig an 
analog system to o sampled data system. If analog break frequencies in the con- 
tinuous system ore greater than half the sampling frequency, — y— , the poles 
corresponding to these frequencies will fold into the primary strip. If a pair of poles 
fold neor the = 0 region it con adversely affect the rigid mode dynamics. 
Choosing a larger samplir^ rate alleviotes this problem, but too large a rate can 
cause the system to be driven by truncation errors. A happy median needs to be 
selected in order to avoid these problems. In Section 2.0 the sampling rates of 
T = 0.04 secorKis and T = 0.08 seconds are used for the synchronous and multirate 
designs. In Section 3.0 the problem of choosing a proper sompl ing rate is treated in 
greater detail . 

2.3 Digitized Continuous Autopilot 


Low Q Sy.^m 

The continuous autopilot G^(s) = 


(10)^ 


can readily be digitized using Tustin's 


(S + 10)^ 2 Z-1 

method with frequency scale pre-warping. This involves replacing S with ^ - y 

2 / "a ^ \ 

and the analog break frequency wlth-^tarJ — ^ j • 


In order to reduce coefficient sensitivity problems 0^(0 h implemented as two 
cascaded first order filters. The general result for one of the cascaded first order 
fi'ters is as follows; 
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where = analog filfer break frequency 

T = sampling rote 
ST 

z = e . 

With sampling periods of T = .04 and .08 seconds and an analog break frequency of 
= 10 rad/sec. this equation becomes 


T = .04 


D.(z)= ^ 

' 1 - .662925 Iz"' 


T = .08 


D,(z) = 


.2971743 (1 + z'*) 
1 - .40565128z‘* 


The corresponding digital representations for G^(s) 
Tustin's method. 


They ore 


T = .04 


D2(z1 = 0.2625 


(1 + z'S 

1 - 0.2504Z* 


18 9 

^ ^ ^ also follow by applying 


T = .08 


DjCz) = 0.4562 


(1 

1 + 0.3056z' 


in order to complete the digital design a digital filter Is added to compensate for 
the phase and gain introduced by the zero order hold. The compensation used in this 
design is the simple lead-lag filter given below. Placing the pole and zero on the 


Gc(z) = 1 .3846 


1 + 0.3 z 
1 + 0.8 z' 
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iwgoHv* real axis inside Hie unU circle allows good zero -order hold compensoflon 
for. small frequencies (o> ■< ) and gain attenuation from the zero"^rder hold for 

large frequencies ( < 6» < ) • The corresponding Nichols charts for the 

low Q system with T — 0.08 and 0.04 seconds are shown in Figures 2.3.1 and 2.3.2. 


Mox Q System 

in this section the effects of digitizing the analog compensation for the max Q model 
ore discussed. For sampling intervals ofT = 0.04 and T = 0.08 seconds the digital 
filters for the continuous autopilot Gc(S) =0.6 f -y 1 and the compensation in 

r 1 LJV*" >-*1-75 T 

the acceleration feedback loop -G q(S) = ^ — J *17 ~ 5 J follows: 


T = 0.04 

• 

^Tustin's j 

Frequency L 

Pre-V/orping 


2 Austin's 
Frequency 
Pre-Warping 


2 Tustin's 
Frequency 
Pre-Warping 

■ .17.5 1 
s+ 17.5J 

T = 0.08 


2 Tustin's 
^Frequency 
Pre- Warping 

[g,«] 

2Tustin's 

^Frequency 

Pre-Worping 

r ^ 1 

[s^lj • 

• 

^Tustin's 

^Fjequency 

Pre-Worping 


[$+ 17.5J 


= 0.6 




0.161379(z+ 1) 


1613792 - 0.838621 




0.020003(z + 1) 


020003Z - 0.979997 


-] 


■ 

r 0.365028(z +1) 

■ 

[ 1.365028Z - 0.634972 


= 0.6 


r 0. 

L 1.33 


3313S9(z + 1) 


331389Z - 0.668611 


p 0,040021(z + 1) 1 

L 1.040021 z - 0.959979 J 


0.842288(z + 1) 


1.842288Z - 0.157712 




-J 
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The Nichol's charh for fhe max Q syst'em with just the analog autopilots digitized are 
shown in Figures 2.3.3 and 2.3.4. Both frequency responses closely resemble the fre- 
quency response of the continuous system shown in Figure 2.1.4. 


The effects of including the digital compensation in the accelerometer feedback loop 
are shown In Figures 2.3.5 through 2.3.10. Figures 2.3.5 and 2.3.6 show the fre- 
quency response of the max Q system with one of the lag filters of Gq(s) digitized. 

The Nichol's chart of the system with ^^'yg^cy ^ ~ s ~ +~ T~l '* Figure 2.3.5 


Pre-Worping 
r* 


r* . " 

and the frequency response with '* 

Pre-Warping L* ^ 


shown in rigure z.o.o. 


In Figures 2.3.7 and 2.3.8 both lag filters in the acceleration compensation are 
digitized. The compensation for the zero-order hold is given below and is the same 
lead-lag filter used with the lift-off flight model. 

^ZOH ” ^ z + 0.3 


Digitizing the compensation in the acceleration feedback loop illustrates an important 
point which should be kept in mind ^en using Tustin's method to synthesize digital 
fiist order lag filters. If the object is to design compensation for a sample data sys- 
tem, using Tustin's method with frequency pre-warping is the same as designing in the 
w-plane. If Tustin's method with frequency pre-warping is used to digitize analog 
compensation, however, the resulting digital filter is used to approximate the analog 
filter's frequency response. Usually this is a good approximation. As shown in 
Figures 2.3.3 through 2.3.6 the rigid mode gain and phase margin for the portially 
digitized max Q system differs slightly from the phase and gain margins of the con- 
tinuous system. When all of the analog compensation is digitized, however, the 
approximation begins to break down. For T = 0,04 seconds the gain margin drops from 
9.14 db to 6,12 db and the phase margin from 51.16® to 42.69®. In additlon,the loop 
at 8 rodlons/sec. present In Figures 2.3.3 through 2.3.6 is missing in Figures 2.3.7 
and 2.3.8. This results from a pair of complex zeros located inside the unit circle 
near 8 radians/sec. moving outside the unit circle when all of the analog compensotion 
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Figure 2.3.3: 

Nichols Chart for Max Q System with Continuous Autopilot 
Digitized using Tustin's Frequency Pre-Warping Method, T = 0. 

Gain Margin = 8*73 db @ 3,74 rad/sec 
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1.3846 ■ ? ^ used for ZOH Compensation 
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Figure 2,3,5: NIchol's Qiart for Max Q System with Continuous Aotopihot 
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Figure 2*3»9: Root Locus For Max Q System with Continuous Autopilot Digitized using 

TustinS Frequency Pre- Warping Method, T = 0.08 
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B digitized. The root locus of the max Q system with just the continuous autopilot 
digitized and the root locus with all of the compensation digitized illustrate this point. 
They are respectively shown in Figures 2.3.9 and 2.3.10 for T = 0.(%. 

2.4 Digital Autopilot Design for the Low Q Model 

As a starting point in the design of a digital autopilot, fiequency response plots are 
run for the low Q model with the analog compensation reploced with a sampler and 
zero-order hold. The block diagram of the system is shown in Figure 2.4.1. The 
uncompensated frequency response plots o- shown in Figures 2.4.2 and 2.4.3 where 
the loop has been opened at K. From these plots it is seen th:at cor.tpensation . V’cS 
'vould phase and gain stabilize the rigid mode and first body U 'ding mode and in 
addition reduce the system gain at frequencies above 11.50 radians/second is required. 


4. 







Figure 2.4.1: Block Diagram of Uncompensated Low Q System 
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Digital System Design for T ~ 0.08 Seconds 

Based upon studying the frequency response In Figure 2.4.3 the rigid mode is phase 
and goin stabilized with the digital equivalent of the following s-plone filter. 

GO) = 0.4764 (^) 


The digitized version of this filter using Tustin's method with frequency scale pre- 
warping is as follows: 


G(z) = 0.6764 


/ 0.2489 Z + 0.2489 \ 
\ 2 - 0.50219 J 


(- 


0.34237 Z + 0.34237 
Z -0.3153 


The Nichol's chart is shown in Figure 2.4.4. 


It was noted in Section 2.2 that a small sampling frequency will cause analog break 
frequencies to fold into the primary strip. This is the case for T =0.08 seconds. 

The body bending poles at ±49.67 rad/sec. in the continuous system fold into the 
primary strip at ±28.87 rad/sec. for the digital system, near the vicinity of the 
second and third body bending poles. This presents an advantage in designing the 
body bending mode compensation. By introducing a notch filter at 27.5 rad/sec., 
additionol gain margin is given to each of the three higher body bending modes as 
well as introducing phase log to the first body bending mode. The phase lag from 
the notch filter is used to help center the first body bending mode about the -1 point 
on the Nichol's chart. Figures 2.4.5 and 2.4.6 show the exploded 
root locus of the rigid mode compensated system and the rigid mode compensated 
system with a digital notch filter Included. The zeros located at z = -0.580 ± j 0.750 
pull the body bending poles at ±27.54 rad/sec away from the unit circle at right 
angles for small values of the open loop gain K. This makes the gain margin of the 
third body bending mode very sensitive to K. The location of the complex poles to 
complete the digital notch filter are placed along the constant frequency lines at 
z = -0.4738 ± j 0.6446. The digital notch filter then becomes 

- , , 0.8459z^ + 0.98 127z + 0.76040 

D (z) = 5 

z^ + 0.947582 + 0.64 

The corresponding Nichol's chart is shown in Figure 2.4.7. 
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Digital System Design for T = 0.04 Seconds 

This section provides a digital autopilot design for the low Q model at T = 0.04 
seconds. In Section 2.5 the design is combined with the autopiiot design at T = 0.08 
seconds to provide a multirate digital autopilot for the low Q model. 


As a starting point in designing a digital autopilot the phase ' . gain margins are 
set for the systems rig'.' mode. This is accomplished by the digital equivalent of 
the following s-plant filter 

G (,) = 0.0279 ( . 

The digitized version of this filter using Tustin's method with frequency scole pre- 
warping is as follows. 


G (z) = 0.8279 


0.13895Z + 0.13895\ /0.19660z + 0.19660 \ 

-T-nrrnm — ) [ — z 0.605791 — ) 


The corresponding Nichol's chart is shown 


in Figure 2.4.8, 


To increase the gain margin of the three higher body bending modes a digital notch 
filter is used. Figures 2.4.9 and 2.4.10 respectively show the root locus of the rigid 
mode compensated system and the rigid and body bending mode compensated system. 

For the rigid mode compensated system the body bending poles at ±27.54 rad/sec 
migrate outside the unit circle for a small value of K, This is evident from both the 
root locus in Figure 2.4.9 and the Nichol's chart in Figure 2.4.8. In Figure 2.4.10 
the body bending poles at ±27.54 rad/sec are pulled into the compensation zeros 
located at z = +0.4000 ± j 0.8500. This keeps the third body bending poles inside 
the unit circl'-. . The location of the complex poles to complete the digital notch 
filter are placed at z = +C.3500 ±j 0.4873. The complete digital notch filter is 
given by D(z) . 

^ , It- 0.800z + 0.8825 

D(z) = (—2 

\ z - 0.700z + 0.3600 
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Figure 2,4.8; Nichol's Chart of Rigid Mode Compensated Lift Off System, T= 0.04 Seconds 

Gain Margins = 8.69 db @ 5.61 rad/sec 
Phase Margins = 33.6" 

First Body Ber.ding Mode Phase Margins * -57* and +120" 

Third Body Bending Mode Gain Margin = fO.52 db 

33 







! VAC SNA'* ' 





«yM OYw: 


ili 



?l JlIN /4 


Figur* 2,4,10: Roo^ Locus of Rigid and Body Bending Mode Compensated Lift Off System 
T = 0,04 Seconds 
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The Nichoi's chart for the rigid and body binding mode compensated system is shown 
in F^ute 2.4.11. It can be compared vv the rigid and body bending mode compen- 
sated system at T = 0.08 seconds which is shown in Figure 2.4,7. 

At this point digital autopilots have been designed which meet or exceed the design 
requirements imposed on the system by the continuous autopilot. In Section 2.5 the 
multirate autopilots For the low Q and max Q systems are presented. 

2.5 Multirate Autopilot 

Two multirate digital autopilot designs are discussed in this section. First the low Q 
digital autopilots at T = 0.04 and T “ 0.08 seconds are used to form a low Q multi- 
rate autopilot. Thb design will be used to outline a multirate design methodology 
in Subsection 2.7. Next the max Q multirate digital autopilot design, using a zero 
order hold os digital compensation in the acceleration feedback loop, is discussed. 

Low Q Multirate Autopilot 

A block diagram of the low Q multirate system is shown in Figure 2.5.1. The rigid 
mode is compensated at the slow sampling rate to keep the system design out of the 
z = +1 region in the z-plane. This helps reduce coefficient accuracy requirements and 
keeps the system from being driven by round-off noise. The body bending filter is 
used for compensating the higher frequencies and is therefore sampled at the higher 
rate. This yields a higher e>^2 which minimizes the amount of rxiise Folding into the 
primary strip. 

Having previously designed rigid and body bending mode filters for the single autopilot 
the same filters are used for the multirate system. In Section 2.4 the rigid mode was 
compensated at T = 0.08 seconds with the digital equivalent of the following s-plane 
filter 

G (s) = K i— ^ where K = 0.6764. 

\s+8/ \ s + 12/ 
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The digitized version of this filter using Tustin's method with frequency scale pre- 
warping is os follows; 


G (*) = K 




2489z + 0.2489 


) (- 


34237z + 0.34237 

7-0V3153 


-) 


In Section 2.4 the 
D(2). 

D(z) = 


body bending mode filter at T = 0.04 seconds was the following 

2 ^ - 0.800z .8825 

2 ^ - 0.700z + .3600 


Combining D(z) and G(z) and using a K of 0.6839 a multirote design is constructed. 

The frequency response and root locus plot for this system are shovm respectively In 
Figures 2.5.2 and 2.5.3. To get slightly better phase stability on the first body 
bending mode the compensation poles in Figure 2.5.3 ore moved to z = -0.18±] 0.3177. 
This corresponds to z - 0.3 ± j 0.519 for T = 0.04. The Nichol's chart for this sys- 
tem with K = 0.81283 is shown in Figure 2.5.4, 


Examining Figures 2.4.11, 2.4.7, 2.1,3 and 2.5.4 reveals the multirate autopilot 
design for the low Q model has no distinct advantage over the single rate design. 

The four frequency responses all meet the same phase and gain margin requirements. 

In addition the rigid mode compensation for both sampling rotes is well removed from 
the z = +1 point. Therefore a multirate system yields no computational advantage 
over a single rate system. 


Maximum Q Multirote Autopilot 

The continuous maximum Q autopilot design os shown in Figure 2.1.4 is used os a 
starting point for the maximum Q multirate design. The continuous system has two lag 
filters in the accelerometer loop -- one breolcing at 1 rad/sec ond the other breaking 
at 17.5 rod/sec. The 1 rad/sec lag filter is being used to reduce the bandwidth of 
this loop while the 17.5 rad/sec filter is used to attenuate body bending effects. 


The max Q multirate design is very similar to the analog design. The 1 rod/sec 
analog filter is left in the loop to reduce the bandwidth of the sampled acceleration 
signal and a sample and hold operating at a sampling rate of T = .08 sec is used in- 
stead of the log filter at 17.5 rod/sec. The sampler and hold contributes a phase lag 
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Figure 2*5«2: Nichol's Chort of the Multjra^e Low Q System 
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Figure 2.5.4: Nichol's Chart of the Fine-Tuned Multirate Low Q System 
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of 40.1* and a gain attenuation of ,74 db at 17.5 rad/sec and the char- 

acteristics are very similar to a log filter ot higher frequencies. This design simplifies 
the requirements for the acceleration loop compensation while not degrading system 
performance. For the multirate design the forward loop analog compensation is replaced 
with a digital autopilot given by 

T = 0.04 

I Tustin's 
Frequency 
Prewarping 

where K = .33113. 

Figure 2.5.5 gives the Nichol's Chart for this multiloop multirate autopilot. The phase 
and gain margins of the rigid mode duplicate those of the continuous system and the 
gain arid phase margins of the body bending modes have been increased. A definite 
improvement has been made with the multirate design. 


r 4 ^ f .0742207 (z + 1)] 

LtiT-JJ ■ I z = . 85135 ? j 


2.6 Pe rformance Results 

The performance results Involve verifying the analysis results via a hybrid computer 
simulation. Open loop frequency response from the hybrid simulation is compared with 
the open loop frequency response generated by the analysis. The hybrid computer is 
used to generate a frequency response for the closed loop system 
, ^ _ G(s) H(s) 

^ ■ 1 + G(S) hIJT * 

The open loop frequency response, G(s) H(s), is derived from the measured closed loop 
response by the following relationships. 


Open Loop Gain 


K 


K = 
o 




K, - 2 cos 0- + 1 

C C V 


Open Loop Phose 
sin 0« 


sin 0,. 


where K = closed loop gain 43 

^ , D 180-18438-1 

0g = closed loop phase 
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Figure ‘2,5,5: Nlchol 's Chart of Max Q Multirate Systenv 
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This assures that all the poles in the simulation ore not in the right half plane or 
on the j 6) axis For the continuous systems and not on or outside the unit circle for 
the digital systems. The systems that were verified are discussed below. 


Continuous Frequency Response for Low Q and Max Q Systems 

The continuous systems were verified by running the frequency response for the original 
analog control system/vehicle combinations. The simulation results are plotted in 
Figure 2.6.1 and 2.6.2. These figures show the analysis results and simulation results 
agree for the continuous system for all frequency values con'ldered 


Digital Filter Programming Form 

Since all the digital filters programmed were either lag filters or body bending filters 
i.e. they involve attenuation, the direct canonic state variable programming form was 
used to realize these filters. A general second order representation for this program- 
ming form is shown in Figure 2.6.3. The relationship between the filter coefficients 
and the scaled coefficients is as follows: 


S = 


SFZ 

0 

JPE7 

1 


’1 


1 


SFE 

o 2 

1 


SFE 

R. = o_ 

® JFT- 


B X 


SFE 

= S?^ ^ ^ gT" 


1 


D' = D X 


E’ = E X 




G^' = Input 1st Summing Point 
(Summing Gain) 

G 2 ' Input 2nd Summing Point 
(Summing Gain) 


A scaled flow diagram is illustrated in Figure 2.6.4. In the results which follow on 
the digital simulations the real axis pole-zero filters are realized as first order filters 
and complex pole-zero filters as second order filters. 
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Digitized Analog System for the Low Q Model, T =0.08 

The next system considered was the digitized analog autopilot for the low Q system. 
The analog' autopilot was digitized using Tustin's method with frequency scale pre- 
warping at T = 0.08 seconds. The comparison between the analysis generoted and 
simulation generated open loop frequency response data is shown in Figure 2.6.5. 

Again the simulation results agree with the anolysis results for all frequencies con- 
sidered . 

Digital Design for the Low Q Model, T =0.08 

This design was for the low Q model with the analog autopilot repi iced with a sampler 
ond zero-order hold. The digital autopilot design appeared in Subsection 2.4. The 
comparison between the analysis and simulation frequency response dato is shewn in 
Figure 2.6.6. The results agree for oil frequencies considered. 

Multlrate Design for the Low Q System 

The comparison between the analysis results and the simulation results Is shown in 
Figure 2.6.7. This plot shows that the analysis and simulation results agree per- 
fectly ot low frequencies and begin to diverge at higher frequencies. This high 
frequency effect was not noticed for the other digital filters because rhe frequency 
response data was not taken at the higher frequencies. This divergence between the 
simulation and the analysis results was expected because of the time delay due to 
A/D setup time, computation time und D/A setup time in the simulaiion. The time 
delay equaled approximately .006 seconds for the XDS 9300 used in the simulation. 

The effect of a time delay is to add phase lag to the frequency response that is 
given as follows: 

0 = -57.3 <uTp 

The triangle points represent the corrected simulator results when the pure time delay 
effects are subtracted out. Therefore, any discrepancies between the results are 
physically explainable. 
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Figure 2*6.5; Digi^i 2 ed Continuous System Nichols Chart for Lift Off Conditions. 

T = 0.08 Seconds 

O Simulation Derived Open Loop Frequency Response Points 
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Figure 2,6.7: Nichol's Chart for Multirate System with Rigid and Body Bending Mode 
Compensation 

(9 - Simulation Derived Open Loop Frequency Response 
^-Adjusted to Remove Computational Delay 
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Muttirote Des ign for the Max Q System 

The comporison between the analysts and simulation results for the max Q multirate 
design is shown in Figure 2.6.8. The results are seen to agree within the accuracy 
of the simulation for the frequencies considered. 

2.7 Multirate Design Methodology 

The following is a general outline to adhere to In designing multirate digital control 
systems. It follows closely the procedure used in designing the low Q multirate 
digital autopilot. 

1. Root Locus of Continuous System - This allows the designer to see the effects 
of different sampling rates on the folding of system poles into the primary strip. 
No matter where the poles of a continuous stable system are located they will 
appear between ± o/2 for the sampled data system. 

2. Z-Plone Root Locus and Frequency Response at Sampling Rates of Interest 
(Synchronous) - By the very nature of including a sample and hold element 
the root locus ond frequency response of the system changes. A root locus 
and/or a frequency response of the basic system is necessary to design any 
type of system compensation. 

3. Determination of High Update Rate and Low Update Rate Compensation 
Requirements - This involves determining the compensation thot is necessary 
at the various sampling rates to meet the systems design requirements. This 
con include not only the various filters that meet the gain and phase require- 
ments of the system but the machine requirements, such as coefficient word 
length, as well. 

4. Multirate F requency Response and Root Locus - From (3) a multirate sample data 
control system is put together. To get a total system description a frequency 
response and root locos plot are made. 

5. Fine Tuning Multirate Design (Based upon (4)) - This can involve merely 
adjusting the system gain or moving the compensation poles and zeros to get 
the final design 


DO 600n 2 145 O = 
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3.0 MACHINE REQUIREMENT TRADES 


The objective of this section Is to point out the various effects software and hardwore 
requirements Impose on control system performance. Computer wordlength, A/D word- 
length, sampling rate, synchronous versus multirate sampling, as well as how the arith- 
metic operations are performed and what state variable form is used, all affect control 
system performance. These considerations are treated In a general manner In Subsection 
3.1. Subsection 3.2 presents a computer wordlength study which specifically answers 
a number of the questions raised by these considerations for the low Q multirate design. 
In Subsection 3.3 an assembly language program for the rigid and body bending mode 
compensated system for T = 0.08 seconds is presented. The puqpose of this subsection 
is to give an indication of the storage and time requirements needed when a digital 
design Is Implemented. 

3.1 Machine Requirements/Synthesis Approach Tradeoffs 

A block diagram of a digital computer and the necessary peripherals to interface this 
system with both analog and digital devices is shown in Figure 3.1,1. In specifying 
the machine requirements it Is necessary to study the resolution ond noise characteristics 
of the various elements In this figure. For example the A/D element will have certain 
requirements and the computer will have memory, word length, and execution time 
requirements. The implication of these machine requirements will be discussed as they 
affect the synthesis approach tradeoffs and vice versa. 



FIGURE 3.1.1; Hardware in a Digital Control System 
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Coefficient' Quantization 

The accuracy that a digital computer can store a coefficient is dependent upon the 
computer intemol word length. The performance of recursive digital filters can be 
adversely affected due to coefficient quantization. This problem has received con- 
siderable attention in the past by Kaiser til and Rader and Gold [21. The sen- 
sitivity of a filter to coefficient quantization has been shown to be greatest when 
the filter poles approach the plus one point in the z-plane. This situation arises 
when the sampling rate is large lelative to the frequency of the digital filter poles. 
Due to the finite accuracy of a digital computer, the performance of a digital filter 
nray become less like its corresponding analog counterpart rather than more, as the 
sampl ing rate is increased . 

Some examples pointing out what is happening will now be reviewed. A simple first 
order example that graphically shows the problem is the following case considered by 
Rader and Gold [21; 

y (nT) = Ky (nT-T) 4 x(nT) (3.1) 

The pole position for this filter is equal to K. When K is realized as K •+• AK the 
error in pole location is AK„ The importance of this error becomes magnified when 
K is near unity. E.g., the maximum gain of this one pole filter is 

(3.2) 

Thus, as K approaches one the sensitivity of this filter to coefficient quantization 
effects becomes large. As things tom out this situation arises frequently in digitizing 
analog filters. 

For example, the z-transform of a first order filter gives the following G(z); 


max 


G(z) 


1 

Tqr 


G(z) = 


T <o.z 
A 



z - e 


(3.3) 
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The Tustin approximation to a First order lag Filter gives 

I + T 0./2 

A 

It can be shown that In the limit as T * 0 the denominators of both these equations 
reduce to the Form 

D(z) = z - 1 + 6 j^T + o(T) (3.5) 

Therefore, as T ' — 0 the poles oF these equations move toward the point z = 1 and 

the Filters become very sensitive to parameter quantization eFFects. 

The Tustin's approximation to a lead lag Filter gives the following equation: 


(3.4) 

) 


G(z) - 


/ ^ “a/* \ 

I'*’ / 


(z ♦ I) 


-( 


G(z) = 






(3.6) 


z - 


1 - 


"T 




1 + 


(Hr 


This equation is seen to have the same type oF denominator as beFore but now the 
numerator has a zero that migrates as a fonction T approaching the plus one point 
as T * 0». In this type oF Filter both the gain and phase become very sensitive 
to coeFFicient quantization. 

Similarly, for second order Filters it has been shown by Rader and Gold [2] that 
the coeFFicient sensitivity problem is compounded. They use the Following second 
order example to demonstrate their point: 

y(nT) = Ky (nT-T) - Ly (nT-2T) + x(nT) (3.7) 
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This filter has complex poles in the z-plane given by 

±j0 

P = re (3.8) 

where 

r= /r 

0 = cos ^ K/2 \/l” 

Explicit expressions con be derived for the pole position errors caused by quantization 


of K and t. if it 

is assumed that the 

errors are small so that Newton's 

approximation 

can be applied. 

t *6 • / 



Ar = 

"Ir ^ 

~w 


A0 = 


“If 

(3.9) 

Using the expressions for r and 0 these equations become 


Ar = 
A0 = 

J al 

2 r 

? tan 0 

‘ ^^2r sin 0 

(3.10) 


The error in the radius r Is similar to the first order case while the error in the 
resonant angle 9 is seen to be very sensitive at small angles. The resonant frequency 
error can also be seen to be inversely proportional to the sampling period T based 
upon the relationship 

A0 = Aw^T (3J1) 


Research by Kaiser (1] has further pointed out that the bound on coefficient accu- 
racy is dependent upon the order n of the filter as well as the sampling rate. Because 
of the multiplications involved between coefficients, it would take twice as many 
digits of accuracy to represent the coefficients of a 2 nth order filter versus an nth 
order filter. From the previous results, it has been demonstrated that increasing the 
sampling rate requires an increase in the accuracy of the digital approximation to the 
filter coefficients. 
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Thus, o potential benefit can be gained in terms of computer internal word length 
requirements by decreasi.ig the sampling rate. This is especially true for filters with 
poles in the vicinity of the plus one point in the z-plone. Root locus methods have 
been applied to study the migration of filter poles and zeros as a function of sampling 
period T to aid in the proper choice of T. Likewise, by realizing the filters as low 
order elements the computer internal word length requirements can again be reduced. 

In general, the transfer functions of a digital controller are realized in the cascade 
forms 


D(z) = 


m 


n 

i = 1 



K. 

I 



+ b.z + c. 
I I 

+ d.z + e. 
I I 


(3.12) 


The real pole and zero filters are realized as first order elements and the filters with 
complex poles and zeros are realized as second order elements. It r also possible to 
realize these filters in the analogous parallel form. The pros and cons of the parallel 
over cascade realization will be discussed in the section dealing with arithmetic noise. 


A/D Requirements 

The requirements on this element are primarily due to resolution and noise level. The 

resolution capability of the device Is proportionol to the number of binary bits used to 

represent the analog signal. The A/D quantization step is generally chosen so that it 

preserves the basic accuracy of the sampled signal The effect of the A/D device on 

filter performance can be statistically analyzed by considering additive noise at the 

input to the filter computation. Assuming each noise sample is uncorrelated and has 

2 

variance Eq /1 2, the steady state variance of the filte - output due to this noise is 
given by [2] : 




0 

IT 


77 


T / (-T-) 


(3.13) 


where Eq = quantization level 

H(z) = filter transfer function. 
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In general the output noise variance is proportional to the gain of the filter and 
inversely proportional to E, the distance of the filter poles from the unit circle^ e.g. 


o 2 
f 



12 E 


(3.14) 


The output noise for certain filters can also be shown to be 
resonant angle of 9 = T of the filter [2] e.g.. 



12 E sin^ 9 


dependent on the 


(3.15) 


From this equation, it can be seen that a low frequency filter (o>^ 0) can have three 

or four more bits of noise standard deviation than a larger resonant angle implementation 
of the same filter. Thus, when trying to maintain a specified signal to noise ratio on 
the output, one might want to think about increasing T to increase the re sonar : angle 
O^T. Root locus techniques can be applied to analyze the filter pole zero migration 
as a function of sampling period 1 . 


Computer Arithmetic Requirements 

The precision of arithmetic operations is related to the machine internal data word 
length and the method utilized by the computer to perform orithmetic calculations. 
The basic structure of a computer data word is shown in Figure 3.1.2. The first bit 
of the data word is the sign bit which is normally 0 for positive numbers and 1 for 
negative numbers. This leaves n-1 other bits to specify the number in binary form. 
When the decimal point for all calculations is assumed to be between the sign and 
the first bit and the remainder of the data word is considered as a single fractional 
number, this is a fixed point data word format. A fixed point data word has the 
following fixed dynamic range: 

1/2"-^ <R< 1-1/2""’ (3.16) 

Fixed point arithmetic utilizes the full wordlength but has the inherent problem of 
scaling quantities to remain within this fixed dynamic range. For example, in a high 
gain filter implementation us^ng fixed point arithmetic, the input signal must be scaled 
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down to prevent overflow. Therefore, the signal to noise performance of this filter is 
adversely affected. To get around this dynamic range problem, floating point arith- 
metic is typically used. In floating point, the data word I: segmentized such that 
part or the word is the so-called normalized fraction of fhe number while the other 
port is the exponent. The normalization process involves shifting the contents of a 
register left until the contents of the registers sign bit and the contents of bit posirion 
one are not equal where vacated bit positions are filled with zeros. Noimalizing causes 
the original fixed point number to be shifted into the most significant bit positions of 
the data word and also maintains a count on the number of left shifts. A floating 
point data word is shown in Figure 3.1.3 where the normalized fraction is stored in 
the fractional part of the word and the number of left shift is stored in the exponential 
section. A general flc. ting point number is thus given by 

N = ±F X 2"^^ (3.17) 

The arithmetic on most flight type computei-s is pcformed on hardware that Is capable 
of performing only fixed point arithmetic. The rea;ons for using these types of machii.es 
are faster computational speed, increased reliability and cheaper cost. Therefore, it 
is necessary to perform some sort of scaling on the filter input, coefficients, and 
sometimes even internal state variables tc not only prevent overflow but to also main- 
tain maxim.um arithmetic precision. The type of scaling used to implement a given 
filter can effect filter performance in terms of signal to noise ratio. This error source 
for uncorrelated input signals can be considered similar to multiplication roundoff noise 
for analysis purposes. [21] 


Sign 

Bit 


0 12 3 4 


i'.-3n-2n-1 n 


• • 


■ V ■« 1 1 1 I I J 

Binary Point (Fixed Point Arithm.-tic) 


JLl 


FIGURF 3.1.2: Computer Data Word Schematic 
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Sign 


Binary 

Point 


n 


i FRACTIONAL * * 

± 

EXPONENT 

1 « ■ » » 


1 1 ■ 1 


J 


n-l n 


FIGURE 3.1,3: Floating Point Data Word Format 


Arithmetic Noise (Uncorrelated Input) 

The primary sources oF computation noise in a digital filter ore due to multiplication 
rouixIofF noise. Similar to the error model for A/D truncation, these error sources 
can be modelled by noise with the following variance: 



(3.18) 


These noise sources are placed in the filter state variable diagram after scaling 
operations and multiplications. For example, consider the state variable diagram for 
the second order difference equation given below (see Figure 3.1.4): 


D(z) = 


-1 -2 
Oq + a^z + 

n -^T' 

1 + b,z + b 2 Z 


(3.19) 


This filter is seen to have 3 multiplications going into the first summer and three 
multiplicatioru going into the second summer. The output noise variance for these 
multiplications is as follows: 


N 

o 





D(z) D(z'^ z“^ 



(3.20) 


This noise variance can be considered as additive to the noise variance of the A/D 
truncation giving a total noise variance for the combination of A/D truncation and 
arithmetic roundoff. 
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FIGURE 3.1.4: Scaled State Variable Diagram for D^(z) 


It is very difficult to give gerteral ground rules for tbe best method of implementing 
digital filters in terms of state variable programming form and scaling technique. Most 
work to date has shown that the parallel state variable form g>vesgood performance 
results for a wide class of filter implementations. But there are even variations in 
the performance of parallel filters depending upon the seal ing method. Thus, the 
question arises as to why the cascade form is considered at all. The primary reason 
for considering it is that the various filter blocks can be changed without affecting all 
the coefficients for the total filter as happens when dealing with parallel filters. Also 
by choosing the proper state variable form, combining the proper poles and zeros in 
the respective blocks and then cascoding the blocks correctly, performance comparable 
to and sometimes better than the parallel form can be derived. It might also be men- 
tioned that although coefficient quantization problems dictate breaking filters into first 
order elements for real poles and zeros and second order elements for complex poles 
and zeros, the computer software can be simplified by utilizing a second order filter 
as the basic building block in either the parallel or cascade form. 

The only general results in the literature are due to Jackson [21] who considers the 
canonic iecond order form (See Figure 2.4) and its transpose and looks at the per- 
formance of this form for both parallel and cascade implementation. He finds that 
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filter performance is dependenf upon the type of scaling used (either max value or 
RMS) and state variable form. He also shows that the implementation scheme is 
dependent upon the measure of signal to noise performance. E.g., if one is trying 
to minimize the output noise variance (total average power) or the maximum value 
of the output noise, different implementation schemes would be required. 

Summary of Machine Requirements/Synthesis Approach Trades 

As pointed out in Section 2.1, the synthesis approach directly affects the machine 
requirements. For example, a high sampling rate synchronously sampled system will 
require much more computer capability than a slow sampling rate version of the same 
filter. In applications with wideband frequency response requirements, the computer 
requirements can often be reduced by going to a multirate realization. The multirote 
implementation could have a fast update rate front end to prevent folding and to meet 
the high frequency requirements, e.g. body bending filters. The low frequency filter 
requirements, e.g. rigid mode shaping filters, could be realized at a lower sampling 
rate which reduced the arithmetic and coefficient quantization problems. Besides 
giving better performance, a multirate filter can olso free up computation time as not 
all the computations are performed every cycle. 

How the filters are split in terms of computation rate is usually obvious from con- 
sidering the respective filter bandwidth requirements. In wideband shaping filter ap- 
plications, it is sometimes necessary to combine any free log filters with the fast up- 
date rate filters as problems can arise meeting shaping requirements at frequencies 
greater than the break frequency of the lag filter due to the notching characteristics 
of digitized lag filters. 

3 .2 Computer Wordlength Study for Low Q Data 

The computer wordlength study was run using rfie low Q multirate system design pre- 
sented in Section 2.0. Four frequency responses were made using wordlengths of 
39, 24, 16 and 8 bits. The scaled direct canonic state variable form was used in 
modeling the rigid and body bending mode filters as it was in the original multirate 
design. A block diagram of this state variable form was shown in Figure 2.6.4. In 
each run the scaled coefficients were calculated, then truncated. The results of each 
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addition and subtraction were truncated os well os the results of each multiplication. 
Before storage, each state variable was truncated. In all four runs the A/D and D/A 
wordlengths were respectively 14 and 15 bits, including the sign bit. 

In Table 3.1 the results of the computer wordlength study are shown for frequencies 
between 0.5 radians per second and 18 radians per secon- . The results show little 
phase and gain differences among the 39, 24 and 16 bit runs. This is not surprising. 
The poles for the rigid mode filter were shown, in Subsection 2.5, to be at z = +.502 
and z =+0.3153 and the two zeros at z = -1 . For the body betiding mode filter the 
poles were placed at z = 0.3 +j 0.519 and the zeros at z = 0.4 tj .85. In each case 
the poles and zeros are well removed from the z = +1 region. This is the region of 
the rigid mode dynamics. The farther removed the compensation poles and zeros are 
from the unit ciicle encompassed by this region the less sensitive the rigid mode dyna- 
mics become to truncation error. Therefore a computer wordlengtfi of 16 bits would 
accurately represent the low Q multirate design. A small computer wordlength could 
also be used to represent the low Q synchronous design. In the synchronous design 
for T = 0.08 seconds the rigid body compensation poles are the same as those for the 
multirate design. The complex poles and zeros fo: this design are respectively at 
z = -0.4738 +j 0.6446 and z = -0.583 +j 0.750. These complex poles and zeros are 
further removed from the z = +1 region than the corresponding poles and zeros in the 
multirate design. For the synchronous design at T = 0.04 seconds the complex poles 
are at z = 0.35 ±j 0.487 and the rigid body compensation poles are at z = +0.7221 
and z = +0.60679. These poles are well removed from “he z = +1 region. Therefore 
a computer with small wordlength could also represent * n Q synchronous design 
without affecting filter sensitivity. 

As stated in Subsection 2.5 the low Q multirate design has no distinct advantage over 
the single rate design. This is not surprising. One advantage that can be gained by 
going to a multirate system is when the design must meet wideband frequency response 
requirements. The fast update rate is used to prevent folding and to meet the high 
frequency requirements. The slow update rate is used for the low frequency require- 
ments and the arithmetic and coefficient quantization problems. In the low Q system 
design there are no arithmetic and coefficient quantization problems with either sampling 
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TABLE 3.1: Effects of Computer Wordlength on the Low Q 

Multirate Frequency Response 




14.1 199” 


3.8 209“ 


9.49 188“ 


-17 247“ 


14.1 199“ 


-3.8 209“ 


-9.47 190° 


-11.6 161° -11.6 162“ 


-11.5 137° -11.54 140° 


-8.75 114° -8.85 116° 


+14.82 318° 16.28 337° 


■16.8 252° 


14.1 199“ 


-3.86 209° 


-9.49 189“ 


-11.4 161.5“ 


-11.54 140° 


-8.86 116° 


13.97 329' 


-16.8 252° 


= 16.0 -27.7 223° -27.21 224.2° -27.23 218.4° 

<u = 18.0 -32.25 187.8° -36.6 189.8° -36.61 184.9° 


11.03 204.9“ 


-4.82 210.2° 
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rafe and all frequency response requirements can be met at T = 0,04 seconds or 
T = 0,08 seconds. There is folding for T = 0.08 seconds. The fourth body bending 
pole does fold into the primary strip but it folds near the second and third body 
be;iding poles. Here it becomes an advantage instead of o disadvantage in designing 
body bendirsg compensation. 

3.3 Machine Requirements 

The machine requirements to implement the rigid and body bending mode filters for 
the synchronously sampled system are discussed in this sccti-: i. Using the Teledyne 
TDY-214 minicomputer and accompanying instruction set, both filters can be realized 
v^ith 47 core locations. A complete list of storage requirements, broken down into 
instruction total, temporary data storage, constants, variable data storage, multipliers, 
adds, loads, stores and transfers are shovm in Figure 3.3.1. 

26 Instructions 

1 Temporary Storage Location 
12 Constants 
8 Variables 

47 Total Storage Requirements 

Average Execution 
Time/Instruction 


12 Additions 

7 1/2 

sec 

12 Loads 

6 2/3 

sec 

11 Store Accumulc‘or 

6 2/3 

sec 

2 Transfers 

6 2/3 

sec 

15 Multiplications 

31 sec 



FIGURE 3.3.1; Machine Requirements for Filter Implementation 


The filter reo'ized with the 47 core locations is the digital autopilot for the rigid 
end body bending modes sampled at T = 0.08 seconds. The filter is shown below 


0.34237 (z-H) Q.248905 (zt-1) 0.57219z^ -h 0.66374z +0.51434 

z -0.315261 z -0.502190“ 1 ^ 0.947582 + 0.64 


The state variables are defined as follows: 

Xl(k+1) = +0.315261 Xl(k) + [0.342370 + 0.315261 (0.34237)1 INPUT(k) 
OUri(k) = Xl(k) + 0.34237 INPUT(k) 

X2(k+1) = 0.502190 X2(k) + [0.248905 + 0.502190 (0.248905)1 OUTl(k) 
OOT2(k) = X2fk) + 0.243905 OUTl(k) 

X3(k+1) = -0.94758 X3(k) + X4(k) + [0.66374 - 0.94753(0.57219)1 OLTr2(k) 
X4(k+1) = -0.64 X4(k) + [0.51434 - 0.64 (0.57219)] OLrT2(k) 

OLTr3(k) = X3(k) + 0.57219 OUT2(k) 


The instruction set for the TOY -21 4 is of the single address type and has the following 
format: 


INSTRUCTION 


ADDRESS 1 

0 12 3 4 

5 6i 

1 7 

Mill 1 

19 


The tog field is two bits in length and indicates which index register, if any, will 
modify the address portion of the instruction. For implementing the rigid and body 
bending mode filters only one of the three index registers Is used as well as just one 
of the two accumulators available with the TDY-214. The assembly language program 
is shown below. 


FILTER COMPUTATION PROGRAM 


(X) 

means 

contents of address X 

Location 

Inst. 

Tag. 

Address 


01145 

LDX 

1 

1200 

Load index register 1 with (1200) 

01146 

LDU 

1 

1205 

Load upper accumulator with (1205 +(index register 1)) 

01147 

MPY 

1 

1202 

Multiply (upper accumulator) by (1202 ’• findex 
register 1)) 

01150 

ADD 

1 

1206 

Add to (upper accumulator) (1206 + (inde- register 1)) 

01151 

‘ITU 

1 

1201 

Store contents of upper occumulator in 1201 + 
(index register 1) 

01152 

LDU 

1 

1215 

Load upper occumulotor with (1215 + (index register 1)| 

01153 

MPY 

1 

1205 

Multiply (upper accumulator) by (1205 +(Index 
register 1)) 

01154 

ADD 

1 

1221 

Add to (upper accumulator) (1221 + (index register 1)) 

01155 

MPY 

1 

1202 

Multiply (upper accumulator) by (1202 +(index 
register 1)) 

01155 

STU 


1177 

Store contents of (upper accumulator) in 1177 


FILTER COMPUTATION PROGRAM (ConHnued) 


Location 


Address 


1214 

1205 

1220 

1201 

1177 

121A 

1210 

1177 

1211 

1210 

XXXX 


Temporary data location 
2 ]q stored here 
Storage location for OUT3 
Storage location for OUT2 
Storage location for OUTl 
Storage location for the input 
0.57219 IQ stored here 
0.248905]q stored here 
0.342370|q stored here 
Storage location for X3(k) 
Storage location for X4{k) 
Storage location for X2(k) 
Storage location for Xl(k) 
-0,94758 IQ stored here 
-0.64|o stored here 
0.502190^0 stored here 
0.315261 ]Q stored here 
0. 66374 ] 0 stored here 
0.51 434 ] 0 stored here 
O.2489O5]0 stored here 
O.34237O]0 stored here 


Load upper accumulator with (1215 + (index 
register 1 )) 

Multiply (upper accumulator) by (1211 + (index 
register 1 )) 

Add to (upper accumulator) (1177 + (index register 1)) 
Store contents of upper accumulator in 1211 + 

(index register 1 ) 

If (index register 1^ / 0, decrement contents of 
index register by one and transfer control to 
location 1146. If (index register 1) = 0 
continue to next instruction 
Load upper O' cumulator with (1214) 

Multiply (upper accumulator) by (1205) 

Add to (upper accumulator) the (1220) 

Multiply (upper accumulator) by (1201) 

Store contents of upper accumulator in 1177 
Load upper accumulator with (1214) 

Multiply (upper accumulator) by (1210) 

Add to (upper accumulator) the (1177) 

Add to (upper accumulator) the (1211) 

Store (upper accumulator) in 1210 
Transfer back to main program 
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4.0 MULTIRATE SAMPLED DATA ANALYSIS METHODS 

In the literature there are two major approaches to analyzing multiloop, multirate 
sampled data systems. The first approach was proposed by G. M. Krone [11]. It 
involves replacing each sampler, T/N, by an equivalent configuration of advance and 
delay elements sampled at a single rate, T. This allows the already developed trans- 
form method of analysis for single rate systems to be applied to multirate systems. 

The second approach is the state transition method introduced by Kalman and Bertram 
[3]. Defining state variables for the continuous and digital elements as well a? for 
the zero order hold elements allows the system to be described by a set of transition 
matrices. The multirate analysis is then carried out by manipulating the various trans- 
ition matrices. The purpose of this section is to investigate these vo ious approaches 
and to choose the most appropriate method for analyzing both synchronous and multi- 
rate sampled data control systems. 

4.1 Krone's Switch Decomposition Method 

In Appendix A two ''•'proaches using Krone's switch decomposition method for analyzing 
multiloop, multirate sampled data control systems are presented in detail. The first 
approach is Krone's original method. The method is relatively straight forward but 
has the disadvantage of involving a substantial amount of work in separating each of 
the branch variables from the governing loop equations. For more compi icated systems 
which utilize a digital computer for analysis, the amount of work required by the com- 
puter to separate the branch variables to form the characteristic equation is greater 
than the work required to do the actual system analysis. The second method discussed 
is the procedure outlined by Boykin and Frazier [12]. It extends Krone's technique 
by introducing a vector notation for the advance and delay elements as well as a set 
of identities needed to convert the resulting vector/matrix equations to scalar equations. 
The procedure simplifies Krone's original technique by allowing the governing loop 
equations to be easily written down by inspection of the signal flow graph. It offers 
little help, however, in separating the branch variables from the governing loop 
equations. 7 q 
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An extension to Krone's switch decomposition method which partially reduces the 
problem of seporating the branch variables from the loop equations is the method 
developed by Thomas C. Coffey and Ivan J. Williams [13], Their procedure is to 
approach the switch decomposition method in the frequency-domain. By incorporating 
appropriate shifts in the s-domain terms and by judiciously choosing the branch variables, 
Coffey and Williams develop an identity which can easily be implemented on a digital 
computer. The disadvantage of Coffey and Williams method is that, like Boykin and 
Frazier's method, considerable care is needed in applying the identity. Therefore 
the analyst needs extensive knowledge of the method or an intricate computer program 
for setting up the analysis. 

Two relatively recent papers byE. I. Jury [14, 15] have extended Coffey and Williams 
switch decomposition method. Jury shows equivalence between the s domain approach 
and the z domain as well as the power of complex integration in formulating the z 
transform. Jury's contribution allows Coffey and Williams method to be worked in the 
z-domain but it too fails to reduce the computer implementation problem. 

All three switch decomposition methods have one common advantage and one common 
disadvantage. From a block diagram or signal flow graph point of view the problem 
is straight forward and easy to understand. The disadvantage is in the degree of 
sophistication needed to implement the various techniques. For this reason, and for 
the relative simplicity in understanding and in implementing the state variable approach. 
Krone's switch decomposition method was dropped as a tool for analyzing multirate 
sample data control systems. The remainder of this section Is devoted to developing 
the state variable approach as an analytic technique for synchronous and multirate 
systems. 

4.2 Stote Variable Approoch to Multirate Sompled Dato System Anolysis 
The application of state variable techniques to multirate systems was first proposed by 
Kalman and Bertram [3]. By using the concept of state and transition matrix, it is 
possible to analyze both synchronous and multirate systems within the same general 
framework. The state approach also overcomes a serious difficulty of transform ir.ethods 
by giving time response information at all times and not just at sampling instants. This 
approach is also easily implemented on a digital computer. 
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A typical linear sampled data control system block diagram is shov/n in Figure 4.2,1. 
From this diagram it can be seen that the sampled system contains the following types 
of dynamic elements: continuous dynamic elements, discrete dynamic elements and 

sample and hold dynamic elements. The dynamics of these various elements can be 
described by state variable relationships. The state of a dynamic element can loosely 
be defined as the minimum amount of information, e.g., an nxl array of numbers, 
which is necessary to calculate the future behavior of the dynamic elemen^ when one 
knows the input/output relationship for the element and the input to the eieinent. 

4.3 States of Linear Sompled Dato Systems 

A linear sampled data system typ :ally contains the following three types of state 
variables: 

1 . Continuous state variables - The continuous dynamic elements are described by 
ordinary linear differential equations with constant coefficients. The state var- 
iables are defined as the outputs of the integrators that appear In the state 
variable (first order) representation of the differential equation. 

2. Discrete state variables - The discrete state variables are governed by linear 
difference equations. These equations are idealized representations of the com- 
putations that are performed in a digital computer. The computer actions necessary 
to implement these equations involve performing arithmetic operations on sampled 
input signals and stored data, storing the dato necessary for the next round of 
computations and placing the output in a register or storage location where a 

d/a device has access to it. The state variables for the discrete dynamic 
elements are the numbers that must be stored in order to perform the next round 
of computations. These states get updated only during a computation and remain 
constant until the next computation involving these states. The state variables 
are defined as the outputs of the delay elements that appear in the state diagram 
representotion of the difference equations. 

3. Sample and Hold State Variables - The sample and hold element represents, in 
an idealized sense, the operation of sampling various types of information and 
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riien applying this information through a smoothing filter fo the input of a con- 
tinuous dynamic element. So basically this element is a discrete-time element 
that gets updated at sampling instants and holds this value until the next sampling 
instant. The state of this element is equal to its output. 


From the state variable descriptions/ it can be seen that the continuous state variables 
transition in a continuous fashion while the discrete time and sample hold variables 
may be thought of as taking place at the sampling instants only. 



FIGURE 4.2,1: A Linear Sampled Data Control System 


The state of the entire system can then be written as the nxl column vector X 
defined as follows: 


X = 


LX^^J 


(4.1) 


/C . 


where X is a y x 1 column vector of continuous states 
8 X 1 column vector of discrete states 


vD . 

X IS a 


,SH 


X isann-(y + S) column vector of sample and hold states. 


4.3 State Transition Equations 

The state transitions for the various dynamic elements will now be described. Before 
discussing these relationships, it is interesting to look at the time sequence of events 
that transpire In a digital cont. I system in an idealized sense. First of all, It is 
assumed that the discrete dynamic elements and sample and hold dynamic elements 


that get updated at a sampling instant transition in the time instant — *- t|^ . 


Likewise the continuous dynamic elements remain unchanged over this instant but 

+ 


traiisition over the finite interval t. 


Figure 4,3,1. 


k 1 


t < t. 


k + 1 


This is shown graphically in 


73 

D180-18438-1 


DO «000 2 ‘49 ORIC. 4/7 1 


jM 





COMPANY 


Discrete Dynamic Elements 

Transition ^ . fi 

Continuous Dynamic Elements 

Z' Transition 




FIGURE 4.3.1: Idealized Time Sequence of Events 


Under tfie assumption that the discrete and sample and hold variables transition 
instantaneously these two state variable forms can be combined and the state of the 
entire system siriplifies to 


Discrete State Transition 

The discrete time states are transitioned by the following relationship: 

+ S d (.^), d r(y i . I,... r.S (4.3) 

j “ ' 

where T is the time it takes the computer to complete the computations. As 
stated previously it is usually assumed that the computations are done instantaneously; 
therefore, the left side of this equation is written as tj^ giving 


,1, ■'ij ■‘j <’k> ^ '•i, n 4 1 '<'k> 
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Sample ond Hold StQ^e Transition 

The sampling process is assumed to toke place instontoneously. The state transition 
equation For the sample and hold element is thus given by 


», ('k > ' 


n 

2 




1 




= y-i 5 + ,,,,n 


(4.5) 


System Tronsition at o Sampling Instant 

Urwler the assumption thai the discrete and sample and .hold elrmetit- oil transition 
instantaneously at a sampling instant, equations (4.4) and (4.5) can be combined Into 
an overall system transition relationship at a sampling jtant as follows; 


x(tj^ ) = B X (y + f r (tj^) 


(4.6) 


where 


B = 


5 — 


I 


I O 

o'.. 


-^n-(S+y)f-^ 
i 

I 


O |0 

I 


d. 


M 


s.. 


O 


f ^ 


nr(6+y) 


y+i 

c 

V+6+1 

- *n J 


n - (5 + y) 

_L 


Continuous System Tronsition 

The continuous dynamic elements are described by the following first order differential 
equation: 

X® + f^m (4.7) 


QO >ut OMie 4/yt 
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The integral of this equation gives the state transition relationship for the continuous 
dynamic e'e.'.icnts and is given below. 

{\ +») = 0*" It) + J 0(t-a) f^m (o) da 

t. 


(4.8) 


where f ~ ^ 1 ^ 


•n order to simplify the expressions that follow, it is convenient to assume the inputs 
to the cont.' joos dynamic elements are from only sample and hold elements. Since 
these inputs are constant between sampling instants, the following closed form ex- 
pression ccn be written for the integral: 


f Xj 0 (t -a) f' do - X. (0 hj (T ) 


(4.9) 


Therefore, the state transition equation (in component form) for the continuous dynamic 
elements is given by 


<i (V+t) = y] 0*^ I| (t) X (t.) + 2 


'j'-k' • ^ h.' ( t) X. (t|^'^ i= 1 , 2 ,... y (4.10) 

i=i I=K+8+1 


0 < T < \ 


System Transition Between Sampling Instants 

TImj total system transition relationship between sampling instonts can be written os 
follows; 


JL 


(4.11) 
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• ^ • «r«M' 


where 


f 


<V ' 


(i^) 


I I 

^ J 'V 

oil j r o 

O |0| i 


-4 




*<0 “ 


- "k> 


<o 


Stobility of Synchronously Sampled Syste ms 

The stability of synchronously sampled systems can be analyzed by using equations (4.6) 
and (4.11). The complete state transition from time n T to time (nxl) T is given by 


x((n+l)T) = ^ (T) (B X (nT))- f r(nT) 

= l(»(T)x(nT) + h(T)r (nT) 


(4.12) 


where 


♦ (T) 

= $(T) B 



1 l>(T) 

= J(T) f 



’’oking the z-transform of this expression and 

rearranging gives 

r 

>1 

Ft V’ 

1 n 

z 1 -HT) J 

z x(0) + 

1 zl “ <P(T) 

The characteristic equation for this system is 

therefore given by 

det 

z J - t(»(T) 

- 0 
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h(T) R(z) (4.13) 


(4.14) 
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These resulh are summarized in the following theorem [31: 

Theorem !. A stationary linear system is stable if all n zeros, z., of the polynomial 
I 2 I - ^ (T) I satisfy 



for i~ I, 2, •••, n* 

These results are c juivalent to the results using z-transform theory. That this is the 
cose will be demonstrated by working an example problem by both z-transform and 
state variable methods. 


Example 4.1 

Consider the sampled data control system shown in Figure 4.3.2. Taking the z-trans- 
form of this system gives the following results: 


G(z) = z 


I - .-T 


z(T - n- e~^) + (1 - T e~^ - e~ ^) 

(z - 1) (z - e"^) 


(4.15) 


The open loop frequency response for this system can be evaluated by plotting the 
system poles and zeros in the z-plane and evaluating the magnitude and phase from 
tnese singularities to points on the unit circle. The root locus for this system can 
likewise be ascertained. 

Example 4.2 

Consider again the system shown in Figure 4.3.2. The stote variable diagram for this 
system :s shown in Figure 4.3.3. From this figure the state equations for the dynamic 
elements are given by 

= -Xj + Xj 

3<2 = Xj (4.16) 

^3 = 0 

end the state t'^rl•ition relationship at a sampling instant is 

Xg (nT ) = r(nT) - kx^ (nT) (4.17) 
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T = I sec. 



Figure 4,3,2: Sampled Data Control System 



Figure 4,3,3: State Variable Diagram for Sampled System 
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The A and B matrices for this system are os follows: 

/-I 1 o\ / 1 0 0 \ 

=(o 0 ll;B=(o 1 0] 

\o 0 0 / 0 oy 


A = 

Evaluating the relationships 

A . J 


(4.18) 


«I(T) = L'' {(il-A)"') 


or fust numerically integrating the equation X ~ AX over a sampling period, the 
transition matrix )3(T) is 


0(T) = 


( 


e I -« 


1 

0 


T - 1 + e"^ \ 

^ ) 


(4.19) 


0 


(4.20) 


The system transition matrix \lt (T) is therefore as Follows: 

_ /e ^ “ lc(T-l + e ^) l-e~^ 

(T) = $ (T) B = I -Tk 1 

\-lc 0 

The characteristic equation then becomes 

|zl - I = z(z^ t (k (T-I + e"^) -] -e"^) z + e“\l-k (1+T) ) + K) (4.21) 

The open loop poles* for this system can be found by evaluating the characteristic 
equation at k = 0 giving 


r (z-i) (z-e ) 


(4.22) 


The open loop zeros for the system can be determined by subtracting the characteristic 
equation with k = 0 from the characteristic equation with k - 1 giving 


z' 7.(T-1 t-e""^) + (1 - T e'^- e‘^)) 

Thus, the system transfer function is given as shown below: 


(4.23) 


G(z) = 


z(z(T-l +e‘^) +(1 -Te"^-e“^)) 
z (z-1) (z-e ^) 


(4.24) 


* These results are bosed upon the assumption that the system transfer function is of 
the form: 

C _ k G(z) _ k N(z > 

R 1 + k G(z) D(z) t kN(z) 
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This transfer function Is seen to be eqjol to the results os derived from z-transform 
theory. The root locus for this system con be found by solving the characteristic 
equation (4.21) at difference values of k. Likewise, the frequency response can also 
be evaluated for this system. 

The presented simple example gives the impression that the state space approach leads 
to a more complicated solution than analysis by the classical methods. This might be 
the case for this simple example except thot a computer implementation of the clas- 
sical approach would require partial fraction expanding the fond ion G(s) and then 
possibly doing a table look up of stored z-tronsforms. It is then necessary to re- 
combine the 2 -transforms to get the transfer function. The state space approach is 
seen to involve matrix manilupatlons which can be readily carried out on digital 
computers. More will be said in the next section concerning the general app! icabil ity 
of state transition approach techniques to the analysis of more complex systems. It 
might also be mentioned that the state transition approach also leads quite readily to 
optimal synthesis and to the concepts of control labil ity and observability for discrete 
and continuous time systems. 

4,4 Extending the S tate Transition Technique to More Complex Syst ems 
Theorem 1 has been shown to be identical to z-transform results that require the system 
closed loop poles to lie within the unit circle in the z-plane for stability. The use 
of 2 -transform techniques has been limited to studying synchronously sampled systems. 
Modifications and extensions have been proposed to the z-transform but these results 
have been shown to be cumbersome to apply and not opplicable to a wide class of 
problems. The state space approach of Kalman and Bertram can be applied to a 
much broader spectrum of sampled data problems. Also, for problems where a stai ionary 
transition matrix can be defined. Theorem 1 is applicable with the appropriately de- 
fined . It has been pointed out in Reference [3] that a stationary transition matrix 
can be obtained whenever the sampling operations repeat In a periodic fashion. Some 
sample cases where stationary transition matrices arise will now be reviewed. 

Cyclic Variable Rate Sampling 

The sampling operation in a cyclic varloble rate system could be assumed to follow a 
pattern as shown >n Figure 4.4,1, The samples are synchronous but the rate is seen to 
vary in a fixed cyclic manner. 


oo 2 • o r. 


7 ' 
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I I I I I It I II I II I 

0 .4 .6 1. 1.41.6 2.0 2.42.6 3.0 3.43.6 4.0 

Figure 4,4,1: The Sampling Process for a Cyclic Variable Rate Sampler 



Figure 4,4,2:Cyclic Sampling Rate System 


B 0(.4) B 0(.2) B 0(.4) 

/ H 

0 ,4 .6 1 . 

Figure 4,4,3: Time Sequence of Tronsitions 


* 
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'he transition matrix for this system will be stationary if a basic sampling period of 
one second is considered. For the simple system shown schematically in Figure 4, 4,2 
the stationary transition matrix would be computed as follows 


First, define the relationships 

i = AX (4.25) 

x(t^^) = 6 x(t^) (4.26) 

Second, solve equation (4.25) for 0(T.) and then use the B matrix to define the 
following stationary transition matrix (see Figure 4.4,3). 

^ir(l) = 0 (.4) B 0 (.2) B 0 (.4) B (4.27) 

The characteristic equation for this system is given by 

det I 2 I - ^ (1) I =0 (4.28) 


and it is thus possible to conduct both frequency response, root locus and time response 
analysis of this type system. 

Multirate Sampled Data System Design 

A multirate sampled data system is assumed to have more than one sampler where the 
sampling process could be as shown in Figure 4.4.4, Multirate systems wiii exhibit 
a stationary transition matrix whenever the samplers act synchronously in a cyclic 
manner. This will be the case whenever the ratio of the various sampling periods 
is rational. For example, for the samplers in Figure 4.4.5 the period T for which all 
the samples act synchronously is T = 1 second. In general, the shortest period for a 
stationary transition matrix is given by 

^ = y ' where p and q are relatively prime integers. 

The sampling patterns have a least common period given by 

T = T.q = T.p 
r j'' 
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DO 6C 


? t 45 O I G 4' t 


T * 1 sec 



FIGURE4.4.4;SAMPLING PAHFRNS FO^ MULTIRATE SYSTEM 



FIGURE4.4.5;TRANSITIONS FOR MULTIRATE SYSTEM 
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When TA is irrational, it is impossible to write a stationary transition matrix. If 

^ J 

p and q are large integers then almost the same situation arises as it takes a very 
large number of intervals to define the overall stationary transition matrix. 

For multirate systems where it is possible to define a stationary transition matrix, 
the analysis proceeds as for the other types of systems analyzed. For the system con- 
sidered in Figure 4.4,4, the operations that must be performed to define the matrix 
(T) are as follows (see Figure 4.4.5): 

^(1) = 0(1/6) 0 ( 1 / 6 ) 02 0 ( 1 / 6 ) 83 0 ( 1 / 6 ) 82 0 ( 1 / 6 ) b, 0(1/6) B (4.28a) 

where the 8 matrices are given by 

B = discrete transition matrix when all samplers are acting synchronously 

= discrete transition matrix when sampler with T = 1/6 is acting 

82 ~ discrete transition matrix when samplers wiHi T = 1/6 and T = 1/3 are acting 

83 = discrete transition matrix when samplers with T = 1/6 and T = 1/2 are acting 

and 0 ( 1 / 6 ) is the continuous transition relationship evaluated for t = 1/6 sec. 

Using the (T) defined by equation (4,28a), system stability can be analyzed by 
applying Theorem 1. Thus, the characteristic equation for this system would be given 
by 

I 2 I - (T) I =0 


Some example solutions using the state transition approach developed in this section 
will now be compared with the solution technique of Krone for the system shown in 
Figure 4.4.6, 

State Transition Solution 


The state relationship for this system between sampling instants isos follows from 
inspecting Figure 4 . 4 . 6 . 




= 0 
= 0 


(4.29) 


85 

D180-18438-1 


00 9000 J ’41 O 5 4' ? » 




TMK 

or in matrix form 
X = Ax 

where 

0 0 10 

. 0 0 0 1 

A = 

0 0 0 0 

0 0 0 0 
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(4.30) 


(4.31) 


The state transition relationship when both samplers act synchronously is given by 


X (nT. ) = B X (n T.) 


(4.32) 


vrfiere 

10 0 0 

0 10 0 

® " 0-100 
10 0 0 


(4.33) 


Likewise, 

described 


the state transition relationship when only the fast sampler is acting is 
by 

X (nT.^) = B, < (nT.) (4.34) 


where 



The two B 


'i 0 0 o' 

0 10 0 

0 0 10 

10 0 0 

J 

matrices are seen to differ only in the third row. 


(4.35) 


Using the relationship 

0 0/N)- l"’ [s I -a]"’ 


T/N. 


(4.36) 


or just numerically integrating equation (4.30) over a period TAJ/ 0 (T/N) is found 
to be as follows: 
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0 (T/N) = 


1 

0 

0 

0 


0 T/N 0 

1 0 T/N 

0 1 0 

0 0 1 


(4.37) 


For this example, assume that N = 2. The stationary transition matrix for this 
sih>ation involves the transition shown in Figure 4.4.7 which is as follows: 


«(T)=0(|) B, 0 (I) 6 


(4.38) 


Performing these matrix multiplications, gives the following <^'(T): 


^(T) = 


2 

T 


-T 
]- 
-1 

-T/2 


0 

0 

0 

0 


0 

0 

0 

0 


(4.39) 


The characteristic equation is given by det j z 1 - i/r j = 0 and using (4.39) this 
relationship is as follows; 

2^ [22+(^-2)z + 1 + ^ T^] =0 (4.40) 


Krone's Time Domain Sampler Decomposition Solution 

It is instructive to rework the example of Figure 4.4.6 using Krone's technique . For 
N = 2, the fast sampler must be replaced with two parallel samplers acting at the 
slower rate but with appropriate time shifts in their channels so that the sum of the 
two is equivalent to the single fast san:pler. The equivalent system is shown in 
Figure 4.4.8. 
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Th« blocks labelled G($) and H(s) reo' 2 senf the rombination of the zero order 
holds and the original transfer functions and ore defined as follows: 


1 -*T 

Gfj) - ' V* - 

s 

H(s) « ' - y -- 


(4.41) 

(4.42) 


From Figure 4,4,8 the system transfer function In the z-domain can be written as 

G(z) 

.TVS , r:-sT/2 --J 


wr " I + G(z) H(z) + Z G($)] Z [e' 


where the characteristic equation is found by equating the denominotor of this 
expression to zero. 


The z-transforms of the vorious terms In equation (4,43) c.? given below: 


«2) - 

z[.-T/2hw] = 

Plugging these expressions into the denominator of (4.43), the churacteristlc equation 
for this system becomes: 

( ^T" ^ J T" =0 (4,44] 


This equation is seen to equal '•^uation (4,40) (the state tronsition derived characteristic 
equation) except for the two extra zeros in (4,40). Tnese extra zeros result from the 
discrete states x- arid and correspond to th 'wo z*::o columns 't: the matrix tA(T). 
Based upon the tact that »/r(T) is beirsg used to define the stationorv t.ansition refotion- 

shfp 

X ( (k + 1) T ) = \^(T) X (kT) + h(T) • (kT) (4.45) 
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The zero columns of ^^r(T) imply thof the states x* and x. need not be known at the 

rti j 

k sampling instant in order to Find the state vector at the (k + 1) sampling instant 

(31. ThuS/ Xg and y , can be dropped from (4.45) and the evaluation of the deter- 
minant for the characteristic equation can be greatly simplified by reducing the order 
of the matrix jz I | . It most be pointed out that ^(T) could not have ex- 
pressed as a simple product of matrices if x^ and x^ had bee-.i dropped before obtaining 
the final result. Therefore, only after ^(T) has been computed can simplifications be 
made. 


4.5 Multirate Sampled Data Analysis Methods 

The following technique' were considered for analyzing multirate sampled data control 
systems: 

1) Time-domain switch decomposition mefood of Krone [11] . 

2) Vector operator signal flow graph technique of Boykin, Frazier, et al [121. 

3) Frequency domain approach of Coffey and Williams [131. 

4) Z-plane technique of Jury [151. 

5) State variable approach of Kalman and Bertram [31. 


The state transition approach of Kalman and Bertram was found to be the method with 
the most general applicability and was also found to be amenable to digital computer 
solution. 


It was demonstrated for systems described by stationary difference equations that the 
characteristic equation is given by the following expression; 

I z I - ^ (T,^) I =0 


where \[i (Tj^) = 

0(T.) = 
®i " 


n 0 (TpSi 

i = 1 

Transition relationship between sampling instants 
Transition relationship at a sampling instant 

I 

£ T. the shortest period fora stationary transition equation, 
i = 1 ' 
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A stationary multirate system is stable if the n zeros of the characteristic equation 
are less than one i.e., j z. { ^ 

If a system is nonstationory, the previous expression no longer holds and one is 
farced to return to the basic definition of stability to analyze these systems. That 
b, a Steele system will have bounded outputs for all bounded inputs. Applying this 
definition of stability to nonstationary sampled data control systems implies that the 
product of the transition matrices approaches zero os the number of products gets 
large. This can be formalized' as follows: 

^(T^) (Tj) ... (T|^) 0 os N » 

From a control system designer's standpoint, this approach to stability is not very 
enlightening as it says nothing of relative stability. 



5.0 MULTIRATE ADAPTIVE CONTROL SYSTEM TECHNIQUES 


This section discusses the results of a literature search, and study of adoptive control 
schemes that appear applicable to space shuttle. Gibson [311 has shown the classical 
adaptive control problem to be separable into (a) parameter identification, (b) state 
estimation, (c) control. The parameter identification task is by far the most difficult 
problem to solve; therefore, it will be emphasized in the work that follows. 


5.1 Online Technique of Lobbia and Saridis 

The first technique to be discussed is due to Lobbia and Saridis [7] . The., consider 
the problem of online identification and control of a stochastic process that can be 
expressed os a linear model with unknown parameter coefficients. This adaptive con- 
trol system minimizes a quadratic performance criteria in an asymptotic sense. Figure 

5.1.1 shows o block diagram of this controller. The system is required to be: (a) linear, 
(b) time invariant, and (c) closed loop stable, and is assumed to be completely con- 
trollable and observable. The disturbances are assumed to be independent and Gaussian 
with zero mean. A known Gaussian sequence with zero mean is added to the control 
sequence to act as a system probe to facilitate identification. 


For the n-dimeusional system, and r-dimensional observation, the input-output in- 
formation is used to define "auxiliary" vectors. These (gr)-dimensional (g < 2n) 
vectors contain information necessary fora new system parameter estimate: 


9 


T 

9.' 

I 


= (0.) contain information on current parameter estimates (gr x gr) 

= fh.T B , h.°^ A B , ..., h° A B 1 

I O I o O I o oJ 


Y. , = [/ (k 1 3 - 1 ), yT (k + 3 - 2 ), .. . y^ (k)] (gr x 1) 

' ^ vector of gr known input noise samples. 


= [o^(k+g-l, u^(k+g-2), ..., u^(k)j (gr x 1) 
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DELAY 


FEEDBACK CONTROLLER 


BoO'l 


MATRIX 

TRANSFORMATION 


IDENTIFIER 


SYSTEM DESCRIPTION: 


S (k+1) = As (k) + B 0 (k) + y (k) + Dv^k) 

(mxl) 

2 (k) = Hs (k) + V (k) 

(mxl) 

o (k) = CONTROLLER 

= PAP"^ B^ = PB, Ho = HP"^ 

(rxl) 

FIGURE 5.1.1: SELF-ORGANIZING CONTROUER IN CANONICAL FORM 
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Assuming that the system is controlled In a stable manner, then the above vectors 
provide a new estimate of the system parameters as Indicated in Figure 5.1.2. This 
Figure also Illustrates the computational sequence of the Identification algorithm. 

The heart of the algorithm Is tlie irerotlon sequence indicated In Figure 5.1.3. 

Possible Problems 

This identification algorithm presents the following questions: 

(1) Would parameter estimation algorithm converge for a non-stationary system, 
i.e. the formulation only provides for the "asymptotic" convergence of a 
quadratic cost function based on a linear, stationary system. 

( 2 ) Would the identification algorithm converge if the true system is nonlinear 
i.e. would the estimation algorithm converge to a finite dimensional linear- 
model equivalent to the nonlinear system. 

(3) If the estimation algorithm is constrained to some n-dimensional linear 
system, and if the actual linear system is of dimention n > n would In this 
case (a) the algorithm attempt to estimate the "dominant" system dynamics 

(b) would the algorithm be stable. 

(4) The impact of practical computer storage and computation time limitations is 
unknown. 

Any real-world problem is to some extent non-stationary and nonlinear. Computer 
limitations dictate limiting system models (and associated calculations) to significantly 
smaller dimensions than the best estimate of the actual system (i.e. space shuttle 
equations). 

Due to the above potential problems, o high technical risk is assigned to the 
implementation of the identification algorithm as given. 

5.2 Maximum Likelihood Method of Stepner and Mehra 

Before a practical autopilot can be designed much preliminary data on aircraft stability 
and control derivatives must be derived. This suggests some sort of off-line identi- 
fication algorithm [16]. A sotisfactory solution of the off-line identification problem 
can significantly simplify the adaptive autopilot problem, by: 
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FIGURE 5.1,2: IDENTIFICATION ALGORITHM 
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(a) establishing relationships between sensed data and key sfobility parameters, 

(b) simplifying the gerteral model by indicating the dominant dynamics. 

For example, off-line analysis [16] may indicate that explicit on-line parameter iden- 
tification is not necessary for a successful autopilot implementation. Hence, we are 
led to a review of the off-line parameter identification algorithms ofStepner and 

Mehra [16]. 

Reference [16] discusses the application of a generalized identificotion method for 
Flight Test Data Analysis. The method is based on the Maximum Likelihood (ML) 
criterion and includes output error and equation error methods as special cases. 

The development of the Generalized Maximum Likelihood (ML) Method for Aircraft 
Parameter Identification has been motivated by several considerations. First, the 
Maximum Likelihood methods are known to provide better estimates than other methods 
[20]. Second, ML methods are more general and can handle both measurement ond 
process noise. In cases where no process noise is present, and the covariance of the 
measurement noise is known, ML methods reduce to the Output Error Method [T , 18] . 
Similarly, in cases where r>o measurement noise is present, the ML Method reduces to 
the Equation Error Method [18, 19]. Third, ML Methods yield realistic values of the 
variances of the parameters. Fourth, ML Methods can be used to estimate the co- 
variances of the noises. This eliminates the problem of specifying the weighting matrix 
or the covariance of the measurement errors. 

Outline of Algorithm 

The maximum likelihood identification method, as implemented by SCI, is indicated 
schematically in Figure 5.2.1. 

It is assumed that the structure of the system is known. The unknown parameters are 
expressed as a vector 0. A maximum likelihood (ML) estimate of 9 is desired. A 
moximum likelihood function is constructed assuming that the constrained model is 
near optimal, so that the error between the model output and actual system output is 
a gaussian zero mean white noise innovation sequence v . 
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FIGURE 5.2.1: IMPLEMENTATION OF MAXIMUM LIKELIHOOD ESTIMATOR 


The ML estimate 0 is obtained by maximizing a conditional probability density function, 
or minimizing N 

J = 4 s 

where B(i) = H. P. + R. 

'll I 

H = model measurement matrix 
P = state error covariance matrix 
R = measurement error covorionce motrix. 



Minimizing Equation 5.1 subject to the covariance model constraints is simplifiev' by 
a steady state B and Kalman Filter assumption [16]. The problem is consequently 
reduced to finding the roots for 


N 

S 



•-«) 


+ 



3 V 0) - 0 

3^0 


(5.2) 


The roots of Equation 5.2 are found by o modified Newton-Raphsan algorithm. 

This modified Newton-Raphsan algorithm b considered the main advance over previous 
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IdenMflcaHon algorifhms. The key feature is the use of the "information matrix" based 
on the innovation sequence sensitivity to parameters: 



(5.3) 


The parameter update is calculated as follows 

-M-' (-1^)' 


(5.4) 


In conclusion, the Maximum Likelihood method developed by Systems Control, Inc. is 
a combination of three steps: [16] Kalman filtering to estimate the states and generate 

a residual or "innovation" sequence; [10) a modified Newton -Raphson algorithm for the 
parameter estimates and [17] an algorithm to estimate the noise statistics (mean and 
variances of the measurement and process noise) 


Results of Identifying Aircraft Stability and Control Derivatives 

Data in three aircraft were analyzed using this algorithm. The three aircraft data 

were from: 

(a) X-22 VTOL computer simulation 

(b) HL-10 lifting body flight data 

(c) M2/F3 lifting body flight data 

In all 23, 20, and 22 parameters, respectively, were identified. Each identification 
task had its unique problemj, requiring non-standard fix up. 






( 


I 



Problems in Identificat ion 

Most of the problems encountered could be related to attempting to identify parameters 
associated with unexcited modes. These problems are called "identifiabil ity" problems. 
For example an inadequate "reference model structure" may force the parameters of the 
model to account for some major unmodeled effect. 
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Errors in insinimentaHon may result in non-physical parameter values. When control 
input is expressible os a linear combination of the A/C response variables, the in- 
formation matrix (M) should be computer leaving out one or more cf the smallest 
eigenvalues. Each eigenvalue which is left out relates to a singular direction in 
parameter space and therefore indicates a combination of parameters which can not 
be uniquely identified. 

Comments 

The Maximum Likelihood identification program described in [16] is considered as the 
"state of the art." The reported performance is considered exceptional. Neither the 
program nor algorithms are available for inspection and verification, i.e., they are 
not government property. Reference [16] reports on the parameter identification results 
and p^erely describes in general the actual identification program. 

5.3 System Identification Techniqi of Lion 
Introduction 

The system parameter identification method of Lion [22] summarized here appears suitable 
for use in adaptive autopilots [26]. The problem considered is the identification of the 
dynamics of a plant where only the input and output can oe observed. Lion's identi- 
fication algorithm does this by continuously adjusting the parameters of a model of the 
plant so as to null some measure of error (« ) between the model and the plant. The 
two most widely analyzed error measures are (a) the response error, and (b) the equation 
error [25] . The response error system has the advantage that only the i’-put and output 
of the system need be measured. The equation error, on the other hand, is an alge- 
braic function of the parameters, and therefore has the advantage of giving a parameter 
adjustment method according to a "steepest descent" law. However, in equation error 
systems, all state variables must be measured or generated (estimated). Lion's "gen- 
eralized equation error" method possesses the advantage of both types; that is (a) only 
input and outr jt need be measured, and (b) it is an algebraic function of the para- 
meters, therefore allowing a true steepest descent path in oarameter space. Global 
asymptotic stability for a large class of inputs is proven. 


riC '>0 0* 2 * 4^ O 4»'« 
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The dcvetopmenf is for e single input/ single output stationory I inear system. The 
system input signal is required to excite all modes associated with the parameters to 
be identified. The method [22] is applicable to nonlinear systems as long as the 
parameters to be Identified enter linearly; ie: parameter "c" In Equation (5.5). 

y + Ojy * AqY + cy^ * u (5.5) 

It is necessary thot the form of the noni inearity be known. 


The usefulness of the identification method [22] in adap ■ outopllot design is its 
rapid convergence rate. If the parameter variotion with time is much slower than 
the system time constants, satisfactory parameter identification is obtainable [22, 24, 26]. 


Summary of Lion's [22] Parameter Identification Algorithm 

Consider the system to be identified represented by the transfer function 


G(s) 


N (s,b) 

oTJSf 


n-1 

E 

1 


b. s 

I 


i-1 


Y 

1 


o.s 


i-1 


y(0 

wr 


(5.6) 


Let "u" and "y” be the plant input and output, respectively. We want a real tir 
parameter estimator such that (a, b) converges rapidly to (a, b), ie; a.-»a., b.— ► 


The parameter estimate error requiring only the system input and output is derived as 
follows: 


N(s,b) U(s) - D(s,a) Y(s) » 0 


(5.7) 


( 






Y(s) » 0 


(5.8) 
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IF w« replace b. and a. bysot.ie esfimators b., a. we have a measure of the error In 
the estimates as follows: 

A A (5.y) 

N (s,b) Y(s) - D(s,a) U(s) Error 

The above N($, b) and D(s, a) are not realistic "state variable filters" since they 
represent differentiation of y and u. Therefore a low-pass filter, H($), is introduced 
as follows: 


H(sJ [N(s,b) Y($) - D(s,a) U(s)] - Err 


(5.10) 


Next polynomials N($, b) and D(s, b) are generalized by replacing by 
"independent" polynomials in s. 


I S'* ■ 

A A ^ A , V ^ 

N(s) = N(p.($), b) = 2 ^ p.(s) t>. 

' i=I 

n-1 

D(s) = D(p.('), S) = Z Pi(*) 


(5.11) 


(5.12) 


The polynomials p.(s) are constrained as follows: 

(a) p.(s) ore "independent", that is there exists no constants c„ such 
that p.(s) = c.. (p.(s) for all i / j . 

(b) H(s) p.(s) are "low -pass" 

With the introduction of the polynomials, pj(s), the correspondence between (a., b.) 
and (a., b.) may no longer be one-to-one. The relationship between (a., b.) and 
(S., ^.) can be derived by comparing coefficients of like powers of "s" in 


D(s) = D(s) and N(s) * N(s) 
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Finally using Equation (5.11) and (5.12) we have the general’zed parameter estimation 

error, £ . 

H(s) [n(s) U(s) -S(s) Y(s)] * C (5.14) 

where H($) p.($) = M. are the "state variable filters." Next define the following 
"state variables." 


y. - H($) p.(s) Y(s) = M. Y($) 

Uj = M. U(s) 


(5.15) 

(5.16) 


Using these definitions Equation (5.14) becomes 
n-1 n-1 

A 

°i ^ 

2 

Having defined the error C a performance criterion P(E) = ^ /2 is selected. The 
method of "steepest descent" or gradient to the P(0) surface leads to the parameter 
estimate odiustment low; 



A 

u. b. = 

I I 


(5.17] 


and 



A 

b. = -k C u. 

I • 



+k e y. 


(5.13] 


( 5 . 19 ; 


The 6=0 surface describes a hyperplane in the space spanned by the parameter 
estimates (a., ^.). The desired estimates of (a., b.) correspond to a unique "match 
point" in the hyperplane. Parameter estimation convergence to this unique point In 
the ( 6 » 0) -hyperplane is guaranteed by excitation of all variables (y., u,). Or In 
other words, the system input must excite all system modes in order to guarantee con- 
vergent estimation of parameters. 
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Extensions: 


Preliminary investigations on the “state varidste Filters" for Lion's algorithm, indicate 
a relationship to the "observer," state estimation filters [27J, This state varioble ob- 
server relotionship allows for rhe exteieion of Lion's algorithm to systems containing 
process ond meosurement noise . Lion's approach to non-l inear systems may provide 
an extension to the "observer identifier" approach of Reference [281. 

Conclusion 

Lion's [22, 231 parameter identification algorithm appears applicable to real time auto- 
pilot applications [261. State estimation filters ("observers") Ref. [27] provide guidance 
in the selection of the "state variable filters" for the identification algorithm. 


The system input must have the necessary harmonic content to excite all system modes 
for the algorithm to converge. This appears not to be a practical problem since ran- 
dom disturbances can provide this excitation. 

5.4 Extensions of Lion's Technique (Luder's Technique) 

Introduction 

Lion's [22] parameter identification technique was recently rediscovered in the 
literature [291. The re-interpretation of Lion's results [29] in the light of observer 
theory [272 provides a minimal stote and parameter estimator. The development is 
for a single-input single output case; extension to the more general is straight 
forward [301. 


Sum mary of State and Porameter Estimator 

Given the system transfer function (5.20) with unknown parameters {o , 0) 


G(s) = 



n 


n 

s 



1 


- 1 



(5.20) 


oe (SC3 1 1«9 OMie. o'ri 
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This transfer function is expressed in terms of known design parameters ( X) gives 
Equation (5.21) 

Jl, »»• 

bi + E ttx” 

G(s)= ? 


(s 


- a,) - E 


a. 

I 


(5.21) 


2 s + X. 


The relationship between {a ,p) ond (a, b) in terms of (X) con be derived by 
equating (5.19) and (5.20). 

The state vector equations (5.21) are motivated by the desire to associate the output 
of Lion's nj state voriabte filter with actual system stotes. See Figure 5.4.1 for the 
block diagram representation of Equation (5.21). 


w « 


. 0 


(5.21o) 




1 ? 


i I 


\ ^ 


% 

/. 

V = 



Vr.F 

• 

T 

n 

II 

o 

II 

•••f b ) 

ft 


V + 


r =(1, 1) 

F — diog ("Xor •••/ “X ) 
* n 


Wi * u, v^ = y 
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The "observer" 1271 equation based on Equation 5.21b with w as input is given by 


r . F 


. 1 J 


b w + K (y - ^ 


if the observer ga!n is as follows: 


c — •••» o) 


(5.23) 


Thus we hove the odoptive observer of Reference tNH* Finally using the "steepest 
descent" RZ parameter adjustment law. 


y. v.t 

I I 


S. w. e 

I I 


(5.24) 


A A 

e=y_y=y_v, 


We have the parameter and state estimation algorithm E9]. This algorithm is 
summarized in Figure 5.4.2 . 


Lion's [22] development is more general than [291. Specifically non-stationary, noisy, 
and nonlirsear problems can be treated by Lion's [221 method but not by the method of 
[291. Analytic extension of [291 and computer simulations are considered to be poten- 
tially fruitful areas of effort [291. 

Example: 

Given: A third order system with unknown parameters (o., b.) 


G(») = 


— 


(s-o,). ^ 


(5.25) 
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Find: The adaptive ccntrol law that fixes the closed-Iuop characteristics, e.g. 

”” suppose it is desired to set the closed loop poles at 


( -r,, -fj, -Tj) 


(5.26) 


Solution: Since the system eigenvalues are invariant under similarity transformation 

■ any convenient coordinate system can be used (see Figure 5,4,1 ), 

= V, = y 


*2 ' V2 * S”: 


(5.27) 


'3 “ Vs * S'“3 


The' dynamics matrix in terms of (z) becomes 

A 
r 




' • 

- z# 


w 



(5.28) 


The adaptive control law is derived by assuming: 


“ ” •') "2 '=3> ( *2 

^3 


(5.29) 


Then Equation (5.28) becomes 


2 = Az + bu = (A + b c^) 2 


(5.30) 


The adaptive gains (c^ , C 2 / ore determined from Equation 5.30 by equating 
coefficients of like powers of $ in Equation (5.31). 

det 1^- I $ + A + b ] = • (* + rp (» + (* •’3) (5,31) 

no 
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I I 

o 
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I 


replacing (a, b) by (a, ^) wii hove the adaptive gain elements (c,/ c«, c,) as a 
function of (a, b, X / r). 

The above algorithm is considered as rtie “state-of-the-art" in on-line adaptive 
control system techniques* 


5,5 Extartsions to Luder's Technique 

An extension to Luder's [28] parameter-state estimation algorithm is presented in this 
section* This extension is considered superior to the Formulation given in [28] 
because of the greater "state variable filter" output separation [22], This greater 
filter separation results in better parameter observability, and hence better algorithm 
performance* 

Brief Development 

It is assumed that the completely observable system can be described by an n-th 
order time -invariant differential equation* For the sake of simplicity, rtie riow can- 
onical form is derived only for the single-input single-output case* Nevertheless the 
extension of this canonical form to the mul ti-input case is straight forward. 


Given: A stable stationary observable system tronsfer function with unknown 
parameters (oi,0) 


G(s) 


1 


(5.3 


s" + 


E 

1 


Ftrxf: Convergent parameter and state estimator 
Solution: 


Restriction of Lion's [22] "state variable filter" to a simple form leads to a state 
estimate (observer) relationship. The transfer function (5.32) can be expressed more 
conveniently in terms of known parameters ( X.) as follows: 


DO cooo ii4S elite* 4/7t 
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■ ■*»r 


G($) 


n-1 

‘‘n ^ E b. 


I ! 

n»r 


(5.33) 


$ + »^ \ ^ V' M. O; 

n ^ I • 

1 


where 


M. = 

I 


n-1 

n 


1 


(s- X.) 


"In expression (5.33) the (a, b) are now the parameters to be identiRed. The trans- 
formation relating (a, b) to ( a , 0 ) involving ( X.) can be derived easily by equating 
coefficients of like powers of s (see appendix B). The form of Equation (5.33) is 
motivated by convergence requirement, and results from the replacement of rfie in- 
tegrators in Figure 5.5.1b by "lossy integrafors" Figure 5.5.2. 


The new canonic form is as follows: 



(5.34) 


(5.35) 
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fS COMPANY 

) 

and where 

J « (b^, ... b„) ; v„ = y = system output 

w = u “ system Input 
n ' 

Figure 5.5.2 summarizes the above state variable form, note that if all \. = 0 then 
(«/^ * (a, b) and we have the "expended controllable canonic form" indicated in 
Figure 5.5.1b. 

Consider Equation (5.35), treat w as a system input, then the state "observer" [27J 
equation is given by Equation (5.36): 
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Next select the "adaptive observer" gain 

- (0, ... 0, 1, 


(5.3T)** 


and substitute (a, b, w) for (ej, b, w) In (5.36) yields the identifier (5.38), which 
is summarized in Figure 5.5.2b. 



The 'kteepest descent" parameter cdiustment !dw [22] adjusts the parameters towards 
the e * (y - ^ = 0 surface. From Equations (5.35) and (5.38) we have; 


nrl 


n-1 


e = -a. V. ^ 53 °j V. > (o^ + X^) y + ? + 0^ y + b^ w - b^ C 


The normal to the e = 0 hyperplane in the space spanned by the parameter estimates 
(a, b) is in the direction of "steepest descent." Define F “ 1/2 (e) as a metric to 
the e = 0 surface, then 


s, . .)•, 


9F 

da, 


y. e V, 


(5.40) 


A 

bj * - 8j 


3F 


at 


- 5. (“ e w.) » + 8j e W| 


(5.41) 


I 

** A minimal stotionory (non-adoptive) observer results if we select Ic = (0,...,0,1,-X ) 


.0 »P0S 0«IC 4'Tl 
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Th« parameter adjustment law completes tfie adaptive observer identification algorittim, 
as su.nmarized in Figure 5.5.3. The above parameter adjustment lows can olso be 
developed os a convenient selection satisfying the estimate convergence proof given 
in the next section. 


Piroof of Estimate Convergence 

The Ljapunov proof developed here is based on [281. This Is considered more elegant 
than the convergence proof given in [223 , since no hueristic argument on hyperplane 
motions is required. Before proceeding with the proof it Is convenient to recall rtie 
superposition theorem; that is the solution to (5.34) in terms of the solution to (5.32) 
Is given by (5.43). 
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CIcarl/ If \ > 0, then the Influence of ony initlol estimate error (Aw(O)^O) will 
decay to zeri os t, — oo . 

Stmtlarly the solution tof 5. 35) in termsofthe vector 2 ~ (Vj, ... y) '** 


V (t) * V (t) + exp (t)j A v(0) 

where v * y, Av(0) = v(0) - y (0) 

n 

Using equations (5.43) and (5.44) in (5.39) results in 

e « - (o^ - a ^ exp (-^^ t) A v(0) 

“Xj,e + (b^ - b^) w_ + b^ exp (~^ t) Aw(0) 


(5.44) 


(5.45) 


iivhcrc o “ (o| f • • • / 



* (|y 

*>n>- 


Now select a Lypunov function 



V= 1/2 e^ + 1/2 IIao 

11/ * 

1/2 |l4b|lj 

n 


.2 

» + 1/2 ^ 

(®i * 

L + 1/2 V ' 

i'i 



I * I 


y A .2 

(^i - '■0 

“T- 


{S,4St 
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Using the stotlonory system assumption = b = 0) we have; 

A 

b. - b. 

I I 

I 1 - 

inserting the "steepest descent" parometer adjustment low, Eq. 5.40 and 5.41, we hove: 

” ” A A 

V = ee + e ^ (o. - a.) v. - e ^ (b. - b.) w. 

1 1 

( 5 . 48 ) 

= ee + e (a^ - a V - e (b - b) w 
Substitute Eq. (5.45) into (5.48) results in 




V = -A. 


^-a^ exp t) Av(0) + b^ exp (.^t) Aw(0)j 


(5.49) 


Recall thot the autonomous or "unforced solution to (5.34) and (5.35) is: 

(>A A^ Yfi) = «xp v(0) + exp (-^0 w(0) (5,50) 


0 os t 


00 if the system to be identified is 


Clearly (y^ + A^ y^) 
stable; i.e. strictly LHP poles. Hence for a stable system (5.49) becomes negative 
definite, therefore parameter and stote estimate convergence is guaranteed. 


implicit in the above is .he requirement thot the input u contain at least n-distinct 
frequencies (real or complex) f22J. 


Convergence Proof for Non-Stotionary Systems 

It is instructive to consider the non-stationary state and parameter estimotion con- 
vergence. These considerations yield guidelines in initial design parameter selection. 
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Consider tfie Hnie derivative of Equation 5.<r6 with a 0 and b \ 0. Carrying 
out steps analogous to those resulting in Equation 5.49 results in Equation 5.51. 

"a ” 

n / a. - S, \ t 

V 

n j 

I ' ' I 




(5.51) 


+ • ^-a^ exp {->^0 Av(0) + exp (-^t) Aw{0) j 


Equotion 5.51 indicates that maximum (X,5,y}are desirable For maximum state and 
porometer estimate convergence rate. A set of large A. > 0 provides for ^ id state 
estimation convergence. A set of large (8,Y) accelerates tne "steepest descent" 
parameter adiustment law. 

A limitation on parameter estimate rate is presented by the frequencies of the system 
to be identified. To see how this comes abour, consider the highest natural frequency 
associated with parameter a.. This natural mode determines the rate of change of 
V. o., which in turn governs the "observability" o' any error in estimate and hence 
removal of the error. 

In conclusion the practical upper bound on (5/X) >0 is dictated by the unknovra 
system natural frequencies. The maximum (A) is dictoted by various numerical con- 
ditioning considerations . 

Conclusion 

It has been shown that for any general input u, the odaptive observer described by 

Equations 5.38, 5.40, 5.41, 5.42 will asymptotically yield the states and parameters 
th 

of on n^ order lineor time-invoriont system [271. This adaptive observer does not 
require auxiliary signals to be fed back into it [271, and hence its implementation is 
very simple (see Figure 5.5.3). 
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5,6 Quost-Adapfive Control Systems 

Mosf of Hie cdoprive control techniques discussed to dote ore still in the research stage 
and present a high technical risk. Therefore, it is necessary to consider proven quasi- 
adaptive systems for a space shuttle type application. The quasi-odaptive scheme to 
be discussed has been applied successfully by Boeing on the U.S. Air Force SRAM 
missile. This is an oerodynomicalty controlled missile that operates over a dynamic 
pressure range of 240 to I inch. The static stability of this missile expressed in terms 
of fins unlocked natural frequency ranges from 14 rad/sec unstable to 10 rad/sec stable 
The structural modes also limited the design of the stabilization compensation for the 
missile . 

The basic autopilot was an attitude type similar to the system on shuttle. The missile 
is guided by an onboard inertial measurement unit so that knowledge of the current 
flight condition is available. The guidance system provides an attitude error signal 
to the flight control system. 

The flight control system design mode use of the aforementioned attitude error signals, 
along with rate gyros in all three axes, forming an essentially conventional ottitude 
control system configuration. Control fins were operated by proportional hydraulic 
position servos. The system was unconventional in that dynamic compensation networks 
were designed with four selectable sets of characteristics to accommodate the wide 
range of flight conditions. Network selection vras determined on the basis of on- 
board computed dynamic pressure, obtained from inertially sensed total velocity and 
altitude. The compensation network design required careful optimization to minimize 
the number of networks required, within the constraints of missile static stabil ity and 
structural bending coupling. This network optimization was accomplished through 
trade studies of gain and phase margins vs. compensation network parameters. 

This design approach is amenable to well established analysis and simulation tech- 
niques, permitting a thorough understanding and evaluation of the system prior to 
commitment to actual flight test. The basic stabilization design was accomplished 
using root locus and Nichols chart linearized stability analysis techniques. Simulations 
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digital/ analog, and hybrid, were used to evaluate nonlinear system characteristics, 
such as aerodynamics and fin actuator behavior. The control system has demonstrated 
unqualified success, performing its intended mission without significant modification 
since the start of flight testing in i969. 

System Description 

The flight control system block diagram for this missile is shown in Figure 5.6.1. The 
missile digital computer is seen to be the broins of this system selecting the rigid mode 
cmd body bending mode filters. The four rigid mode filleti are chosen based upon 
dynamic pressure while the two body bending filters ore selected based upon moment 
of inertia data. 

Flight control gains were computed onboard based on inertiolly measured information, 
along with programmed representations of conh-ol fin effectiveness and missile moment 
of inertia. The control gain computation and network selection logic are shown in 
Figure 5.6.2. 

The onboard determination of flight control gains and selection of compensation 
networks os functions of missile -measured flight condition constituted a quasi -adaptive 
mode of operation. The system adjusted itself according to measured flight condition, 
but not through o dynamic mode of closed-loop sensing in the 'suol sense associated 
with adaptive control. This is a significant advantage when there is a severely con- 
straining closeness between the frequency regime of rigid body control and that of 
structural body bending, as was the case with SRAM. 

The high degree of success of the design approach used on SRAM merits consideration 
of this concept on other appi icotions. It is especially amenable to any system that 
requires an onboard digital computer to perform the guidance and navigation tasks. 

The computer's functions con then be expanded to provide computational services for 
the flight control stabilization task, using the basic information that is generated os 
port of the inertial guidance function. The use of variable compensation characteristics 
and flight control gains based on inertiolly sensed flight condition is not restricted to 
an attitude controller; it is conceptually applicable to normal acceleration control 
systems, or to any technique used to augment stability or derive stability fora 

statically unstable or marginally stable vehicle. 
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This quosi-adaptive design approach is applicable to any vehicle with widely varying 
flight regimes. Therefore, a design similar to this could be applied to the shuttle 
vehicle. 

While *n the mixed reaction-]et/aerodynamic and totally aerodynamic modes the filters 
could be switched as a function of dynamic pressure, the loop gains could be changed 
os a function of fin effectiveness and the mixing of control modes could be performed 
as a function of dynamic pressure. 

5.7 Digitizing Adaptive Control Systems 

The previous discussions have not considered methods of choosing sampling rates for the 

adaptive control system. Two methods of converting parameter adoptive control systems 

into multirate designs will now be discussed. The first involves using the identified 

parameters to determine the dominant response modes of the system. It is precisely 

these modes which dictate the sampling rate requirements from a frequency folding 

effects standpoint. Since the folding effects usually place a higher requirement on 

sampling rate than stability considerations, this technique could be used if analysis 

verified that this was in fact the case. The second involves adapting the sampling 

rate to meet the time varying control system requirements. Assuming that the control 

system was functionally dependent on time varying parameters, a digital implementation 

could require different time varying sampling rates on the various control loops. From 

strictly a digital realization standpoint the sampling rate could become a function of 

the time varying filter break frequencies e .g. a> > N where ^ Ihe 

s — max max 

highest filter break frequency. Care would have to be exercised in blindly applying 
this technique especially if the high frequency behavior of signals being sampled was 
unknown. This situation could be alleviated by analog prefilters on specific signals. 

A possible implementation might involve using both of these methods. 
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APPENDIX A 

THE SWITCH DECOMPOSITION METHOD OF ANALYZING 
MULTIRATE SAMPLED DATA CONTROL SYSTEMS 

U is illushxiMve to work a sample problem using both the original switch decomposition 
method of Krone [11]and the updated version using the vector operation signal flow 
graph technique. First let us find the system transfer runction shev/n in Figure A-1 
using Krone's original method. 



FIGURE A-1 

KRANC'S ORIGINAL SWITCH DECOMPOSITION METHOD 


The system in Figure A-1 is redrawn replacing the T/2 sampler by the equivalent 
parallel combination of samplers sampling at rate T. The equivalent system is shown 
in Figure A -2, 



FIGURE A-2 

EQUIVALENT SYSTEM IN BLOCK DIAGRAM FORM 
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The loop equations necessary to find thv transfer function are written in terms of the 
unknowns B^($) and 82(5) and the output quantity C(s). They are 



(A-1) 

(A-2) 

B2(s)(A-3) 


The Laplace transform of each sampler output is obtained from (A-1), (A-2) and (A-3). 


They ore 




C*^ (A-4) 



Rearranging and toking the z-transform yields 

2(1 + 1/s) B^ (z) + Z^^2 ^ Z (l/s) C(z) = Z(l/s) R(z) 


(A-7) 


z Z,^2 + Z (1 + 1/s) 82(2) + z Zj/2 ^ ^1/2 

* ^'/2 ^ * TSTTf) 


= Z 


(n^ny) 


R(z) 


(A-9) 


where 


= e*^, Z(f(s)) = F(z) and Zy^ 2 ^^i)) = F(z, 1/2) 


(A-10) 
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Applying Cramer's rule to Equations (7), (8), and (9) yields the transfer function 
as a ratio of two determinants. 


C(z) 

wr 


C(z) 

Wf 


Z(1 + 1/s) 

z Z^/2 


Z (1 + 1/s) 
z Zi /2 (1/s) 


Z,/ 2 (l/s) 

Z (1 + l/i) 


Z,/j(l/s) 

z (1 + lA) 


(nrny) ^\/2 (Wtti) 


z (lA) 

z Z I ^2 ( ^ /®) 


(nrrn) ^1/2 (nr^) ^ (irr^) 


z(iA) 

Z,/2(Vs) 

^ (’ ^ TTI^)) 


-(A-11) 


Expend. ng each determinate along the first column allows C(z)/R(z) to be written as 
a ratio of two polynomials. 


C(z)_ 


(t^) * ^ ^(rhn)} 


Z,/2('A)Z,/2 


1 


TTTin) 


) 


1+2Z(1/s)+Z(1/$)Z(1/s)-zZ,/2(1/s' Z ,/2 (Vs) + Z ( 7 ^)) ^ ^ 


Z ^ (ttAt))- ^1/2 h /2 (17^) 

Each of the terrm in the numerator and denominator of C(zV^R(z) can be converted 
to (unctions of 2 by consulting a table of z and modified z transform poirs. 

Vector Operation Signal Flow Graph Technique 

Using the vector operation signal flow graph technique the flow graph for the system 
shown in Figure A-1 can be drawn. The flow graph is presented in Figure A-3. 
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FIGURE A-3 

EQUIVALErOT SYSTEM IN SIGNAL FLOW FORM 
The equations for each output node iirtned lately follow froiT Figure A-3. 


B = (1/s) (R*^ - E^“ - C*^) 


(A-13) 


c - (1/s) 


(A-14) 


where E is the column vector of advances with components E. = e , 

“• N- ' 

i = 1, 2, N and E is the row vector of delays with components 




The Laplace transform of each output node, after sampling, is given by Equations (A-15) 
and(A-16). 

bJ = (1/s) 1*^ (R*"^ - eJ^ " B*"^ - C*^) (A-15) 

= |(l/s) (l/$+l)j* (R*^ - - C*^) (A-16) 


If we shorten our notation by defining 

D*T ^ I |n X N * OA) e"*- } a*' t j E^* (,A) j 

|(1A) (1A + I) ond F*’ 4 |(1A) (lA + >){* 
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where the prime (A-1) indicates the transpose, (A-15) and (A-l$ combine to form the 
fbllowing matrix input-output transfer function. 


C(z) 


F(z) - b^z) d \ z ) g (z) 

1 + F(z) - b' (z) D ^ (z) a (z) 


(A-19) 


Using the identities developed by Soyhin and Frazier 02] C(z,/R(z) can easily be 
conx'erted to a ratio of polynomials. The identities, which are used to break *wn 
vector/matrix relationships of the ferm G, G , and G E^~ intc scalar 
equations are listed in Equations ^-2C| through (\-23). 


(gN^ Y E^“)(z) = (fpV 

' ' P = 0 


'.p>^ <%v 


(A-20) 




(A-21) 


(ye'"')(z) = 

\ - / P= 0 


vVn , > . 

Y (z w ) V 

P P -P 


(A-22) 




M 1 y^/N /, V p T/N, . 
1 X (^p Wp) R (Zp Wp) 


p p 


, _ i2ir./N sT sT/N 

where w-e ,z = e,z.,=e 

p N 


(A -23) 


X, Y and R are scalar functions in z w 

P P 


e 

-P 


^IP 


■ 1 ' 


1 

^p 

V = 

S'J 

• 

-p 

IP 

e 

N-I 

L IP J 


e 

• 

IP 

_ 1/N i2irP/N . 

= z e P 

= 0, 1 
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vT/N , x_ 

X pto ^ ^ 


N - 1 


X^ <Z}s|*^' ^ ^ '*'k^ ^ ** (A -24b) 


To expand (A-19) into a ratio of polynomials let G - 1/s and H - 1/s+l and note 


N-1 

2 e„ V _ - 

_p _p - I 

P = 0 


i N X N 


(A-25) 


Expanding D (z) 


D‘Vz) = 


[! .(■“••!"-) MJ 




v] 


’n-1 / 

= 2 (® i' ) ( 

p= 0 "P P \ 1 t- 


^T/N , X 
G (z w ) 
P P 


(A-26) 


Using Equations (A-21)/ (fi^-22) and (A-2$ 


b' D"’ a = ^ 2 G (z 




/ N - 

(l/N 2 

\ m = 






N -1 G (Zp w ) g"^*^ {z w ) 

= l/N 2 U ^ 


P= 0 


1 + G {z^ Wp) 

132 

D180-18438-1 


4 ' 7 ' 


COMMA'S > 


Finding ~ G (zj^j) ‘ lefting N = P in Equation (A-24b) 

j. N - 1 G H^^(z w ' + G^^(z w ) G ) 

(*M ) = 2 ^ ^A-28) 

N P=0 (^nV 


Subtracting expression (fi>r~27) from (^-28) yields the numerator of the Input-output 


transier hinction 


F^ (zj^N) - b ' D 'a (Zfg*^) = VN 2 


N - 1 


P = 0 


G H (z w ) 
P P 


’ " G (-nV 


Therefore 


C(z) _ 

wr~ 


1 V ' ° V 

N pf-Q 

. N - I G (z w ) 

2 t/n^ ^ 

>-0 '*<5^'^ (-pV 


The expression is valid for any positive ‘nteger N. For our example N - 2 and 
Equation (A-30) reduces to 


C(z) ^ 


Z (GH) + Z(G) Z(GH) - 2 Zi/2 (G) Z]/2 (GH) 

] + 2Z(G) + Z{G) Z(G) - z Z^^2 ^1/2 ^ 


- z Z ,/2 (G) Z ,/2 (GH) 


If l/$ is substituted for G and 1/s + 1 for H Equation (A-31) is equal to Equation 
^-135, the answer which was obtained using Krone's oiTginal switch decomposition 
method. 
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Both Krone's switch decomposition method and the vector operation signol flow graph 
technique hove inherent advantages and disadvantages. Kranc's original method 
is straight forward but becomes cumbersome as the complexity of die system increases. 
The vector operation signal flow graph technique has the advantage of all. wing the 
governing loop equations to be easily vrritten down by inspection of the flow graph. 
The resulting vector/matrix equations however con become hard to deal with even 
though the identities and relationships introduced by Boykin and Frazier [121 makes 
this task somewhat easier. 
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ft .® * 

ow- _i ^ 1 ! 

* + 0^s+03*+a2$+0| 


(B-2) 


Multiplying both numerator and denominator of Equation (B-1) by (s + A.^) (s+AjKs + X,) 
results in a ratio of polynomials of the form of Equation (B-Z) . Comparing coefRcients 
of like powers of "s" results in the following tronsformations. 



Note thot the triangular form simplifies inversion. 
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